Describing systems of interacting fermions by boson models: 
exact mapping in arbitrary dimension and applications 
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We develop a new method that allows us to map models of interacting fermions onto bosonic 
models describing collective excitations in an arbitrary dimension. This mapping becomes exact 
in the thermodynamic limit in the presence of a bath. The boson models can be written either 
in the form of a model of non-interacting bosons in a fluctuating auxiliary field or in the form of 
a superfield theory of interacting bosons. We show how one can study the latter version using 
perturbation theory. Using the developed diagrammatic technique we compared the first two orders 
of perturbation theory with the corresponding results for the original fermion model and found 
a perfect agreement. As concerns the former representation, we suggest a scheme that may be 
suitable for Monte Carlo simulations and demonstrate that it is free of the fermionic sign problem. 
We discuss in details the properties of the bosonic representation and argue that there should not 
be any obstacles preventing from an efficient computation. 

PACS numbers: 71.10.Ay, 71.10.Pm, 75.40.Cx 
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I. INTRODUCTION 

Many body electron systems with interaction tradi- 
tionally attract a strong attention. Numerous interest- 
ing phenomena of the condensed matter physics originate 
from the electron-electron, electron-phonon and other 
types of the interaction and do not exist in the ideal 
Fermi gas of the electrons. 

Generally, models with interactions cannot be solved 
exactly in the space of dimensionality d > 1 and one has 
to use approximate schemes of calculation. The most de- 
veloped method is the perturbation theory with respect 
to the interaction. Within this theory one starts with 
the ideal Fermi gas and calculates corrections to physical 
quantities assuming that the interaction is weak. Some- 
times it is sufficient to compute only the first several or- 
ders of this perturbation theory, but very often in order 
to capture important physics, one has to sum certain se- 
ries. 

The formalism of the perturbation theory for the in- 
teracting electron gas is well developed (see, e.g., Ref. HI). 
Each term of the perturbation theory is represented di- 
agrammatically, which enables one to carry out quite 
complicated calculations. Many interesting phenomena 
like, e.g., superconductivity are not seen in the lowest 
orders of the diagrammatic expansion and show up after 
summation of certain ladder diagrams. Consideration of 
chain diagrams is necessary in the microscopic theory of 
the Landau Fermi liquid and there are many examples of 
this type. 

At the same time, such sequences of diagrams corre- 
spond to physically well-defined collective bosonic exci- 
tations. In the Fermi liquid theory, the chain diagrams 
describe particle density and spin excitations that can be 
spoken of as bosonic quasiparticles. From this and many 
other examples, where a well-defined physical quasiparti- 
cle is represented by an infinite sequence of diagrams, one 



can conclude that the conventional diagrammatic tech- 
nique is not always the most convenient technique for 
the description of interacting many body systems. 

The inconveniences of the conventional perturbative 
methods are especially evident in situations when the in- 
teraction between the quasiparticles is important. This 
happens in one-dimensional systems, d = 1, but mod- 
els in higher dimensions, d > I, used for describing high 
Tc superconducting cuprates^, heavy fermions^, etc., can 
hardly be successfully studied by conventional diagram- 
matic approaches either. These difficulties call for con- 
structing different approaches dealing with collective ex- 
citations rather than with single particles. 

Such an approach called bosonization is well known for 
one dimensional (II?) systems. A huge number of pub- 
lications are devoted to it (see, for a review, Refs. Ijja). 
Attempts to bosonize higher dimensional fermionic sys- 
tems have been undertaken in the past starting from the 
work by Luther— where the Fermi surface of a special 
form (square or cubic) was considered. This idea was fur- 
ther developed by Haldane^ who suggested to patch the 
Fermi surface by a large number of segments and contin- 
ued in a number of publications^^—. All these bosoniza- 
tion schemes are based on the assumption that only small 
momenta arc transferred by the electron-electron inter- 
action such that, after any scattering event, the electron 
remains in the same segment of the Fermi surface. This 
case is, of course, relevant to a long range interaction 
but a generalization to an arbitrary interaction and arbi- 
trary momentum transfer is hardly possible within this 
scheme. At the end one can reproduce results of the 
random phase approximation (RPA) that can be inter- 
preted in terms of quasiparticles but the interaction be- 
tween these quasiparticles cannot be taken into account 
properly. 

Another method based on quasiclassical equations has 
been developed recently for an arbitrary interactioni^. 
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Using this approach anomalous contributions to the spe- 
cific heat and magnetic spin susceptibihtj^ have been 
calculated in all dimensions and new logarithmic contri- 
butions were found. In ID, the results obtained are in 
agreement with those known from renormalization group 
considerations^^ and exact solutions for spin chaina^^ii^. 
They have also been confirmed in 11? by the conventional 
diagrammatic expansion^S where the first non-vanishing 
logarithmic contribution was calculated. However, the 
same diagrammatic method did not exactly reproducG^i 
the logarithmic contributions to 6C/T'^ for d > 1, where 
6C is the anomalous correction to the specific heat. This 
discrepancy has been attributed to the observation that 
not all the effects of Fermi surface curvature have been 
taken into account in Ref. [T^. 

One can conclude from this discussion that the prob- 
lem of the bosonization in c? > 1 has not been completely 
solved yet but such a method is very desirable because 
it would be a new tool for analytical calculations. Re- 
markably, fcrmionic models are shown to be problematic 
for numerical computations as well. The powerful Monte 
Carlo (MC) method suffers from the famous fermionic 
sign probleni22ri2&. The MC method implies that the par- 
tition function can be represented as a sum of terms with 
positive probabilities and in this case, the computation 
time is proportional to a power of the size of the sys- 
tem. This makes the MC method very robust compared 
to other approaches. However, the fermion determinant 
is negative for some paths and, actually, its average sign 
becomes zero in the sampling process. As a result, the 
computation time grows exponentially with the size of 
the system or inverse temperature and the advantages of 
MC get lost. 

We are not able to list here all publications where the 
sign problem was discussed but (except special cases like 
electron systems with attraction, systems with repulsion 
but with half filling and some others) the solution has 
not been found. Moreover, according to Ref. [2^, the sign 
problem should be NP-hard^, which means that its reso- 
lution is almost impossible. We are not aware about any 
mathematically rigorous proof of the NP-hardness of the 
fermionic sign problem but it is anyway clear that it is 
the main problem for MC simulations. 

This problem does not arise in (non-frustrated) bosonic 
systems, though. Therefore, it is quite natural to think 
about a possibility of circumventing the sign problem by 
mapping the initial fermionic model onto a bosonic one. 
Of course, it is impossible to convert real fermions into 
real bosons but one can think of recasting the interacting 
part of the fermionic action in the language of collective 
bosonic excitations. This leads us again to the idea of 
the bosonization. 

In this paper, we suggest a new bosonization scheme 
that allows one to map in the fermionic model to a model 
describing bosonic excitations exactly in the sense that 
the mapping works in any dimension at any tempera- 
ture and for arbitrary interactions. It can be written 
in a form of a model of non-interacting bosons in an 



effective Hubbard-Stratonovich (HS) field. The trans- 
formation from the fermionic system to the bosonic one 
corresponds to using in quantum mechanics the density 
matrix instead of the wave functions. We argue that the 
representation of the model in terms of bosons in the ex- 
ternal field should be free of the negative sign problem 
and can be used for MC simulations. 

It should be emphasized, though, that for finite sys- 
tems used in the MC simulations the new bosonized 
model is not exactly equivalent to the initial fermionic 
one. In order to derive the boson model, we need a 
certain regularization and a bath attached to the sys- 
tem. This is necessary to avoid some uncertainties in the 
bosonized expressions. Only after the regularization is 
carried out we obtain a bosonic model free of the sign 
problem. Completely exact mapping of the fermionic 
model onto the bosonic one is impossible because, tech- 
nically, several electron Green functions with coincid- 
ing Matsubara frequencies and momenta cannot be con- 
verted into bosonic excitations. If we did only identi- 
cal transformations we would not be able to get rid of 
the sign problem. At the same time, one can see from 
the derivation that considering a sufficiently large system 
with the bath any desired precision can be achieved. 

Following an alternative route one can average over 
the HS field before making any approximation. This 
goal is achieved using the famous BRST (Becci-Rouet- 
Stora-Tyutin) transformation based on the introduction 
of superfields^^. This allows us to represent the partition 
function in such a form that the average over the HS 
field can be done immediately. After the averaging one 
comes to a model of bosons with quadratic, cubic and 
quartic interactions. We show how this model can be 
studied with the help of perturbation theory in the inter- 
action terms and compare the first orders with the corre- 
sponding terms of the conventional perturbation theory 
demonstrating a full agreement. 

The paper is organized as follows: 

In Section [Hi we formulate the model of interacting 
fermions, decouple the interaction with the help of the 
HS transformation and map the model onto the model of 
bosons in the external field. We introduce a regulariza- 
tion that allows us to avoid singularities in the equation 
of motion. 

In Section |Ih1 we introduce superfields, average over 
the HS field and derive the model of interacting bosons. 
After that we develop a diagrammatic technique for the 
superfield theory describing the bosonic excitations and 
demonstrate in the first orders in the interaction its 
agreement with the conventional perturbation theory. 

Section HVl is devoted to the discussion of the bosoniza- 
tion from the point of view of MC simulations. We dis- 
cuss in detail how the negative sign contributions origi- 
nate in the fermionic formulation and why they do not 
exist in the bosonic regularized model. Properties of the 
partition functions are considered in the complex plane 
of the electron-electron interaction, which helps to un- 
derstand in details what happens in the procedure of the 
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bosonization. We derive formulas directly suitable for 
numerical simulations and give an example of numerical 
convergence in a case of static fields. 

We discuss the obtained results in Section IVl 
The main results of this paper have already been pre- 
sented in a shorter publication^^. 



II. FERMIONIC MODEL AND BOSONIZATION 

A. Reduction to fermions in the auxiliary field 

Our method of the bosonization is quite general and is 
applicable to a broad class of models. Therefore we start 
our considerations with a general model of interacting 
fermions on a lattice of an arbitrary dimension. We write 
the Hamiltonian H of the system as follows 



H — Hq + Hint, 



(2.1) 



where is the bare part describing non-interacting elec- 
trons, 



(2.2) 



r.r' .(7.(7' 



with a chemical potential ^, a hopping amplitude t^-y 
between the sites r and r' and the spin cr = +1, — 1. 

The term Hint describes the electron-electron interac- 
tion, 



Hint — „ ^ ^ V^.r'C^o-C^'o-''''''^'' 



(2.3) 



r,r' .(7,(7' 



The partition function Z can obtained from H the 
Hamiltonian H using the standard relation 



(2.4) 



Z = Tr exp , /? = l/T, 



where T is temperature. 

In order to proceed with the mapping of the fermion 
model onto a boson on, we will decouple the interaction 
term Hint by integration over an auxiliary HS field. In 
principle, the decoupling can be performed in different 
ways and we should choose one of them. The choice 
can be different depending on the problem studied. In 
this work we are mostly interested in systems with an 
on-site repulsion. We also keep in mind a possibility of 
using the method for MC computations. The latter desire 
enforces us to use real HS fields 0^ depending only on one 
coordinate r. 

At the same time, it would be advantageous to keep an 
arbitrary sign of the off-diagonal interaction Vr^r' ■ The 
on-site attraction is less interesting for us because MC 
schemes can be free of sign in this case^^ anyway. The 
bosonization procedure can be carried out also in this 
case with minimal changes but we do not do this here. 



In order to carry out the HS transformation with the 
real fields (j)r we split the matrix Vry into its diagonal 
Vr^r and off-diagonal Vr^r' parts, 



Vr 



Vr 



Vr 



0. 



(2.5) 



Assuming that Vr.r' is bounded from above, we can 
add a sufficiently large constant V to —Vr^r', ensuring 
the matrix 



(1) 



-Vr^r' + ySr,r 



(2.6) 



to be positive definite. (Its Fourier transform Vq^^^ should 
be positive for all momenta q.) Introducing the coupling 
constant 



Vq = Vrr + V 



(2.7) 



assuming it is site independent, and using the fermionic 
commutation relations we rewrite the interaction term 
Hint as 



n-int ~ J^i„t + -n. 



H 



(0) 



ijit ' int ' 



(2., 



H 



(1) _ 



- y v^'] 

9 / ^ r,r' 



(1).+ 



while the chemical potential /i should be replaced by /i — 
M - (^^0 + V) /2. 

Due a somewhat arbitrary choice of the constant V 
the coupling constants Vq and V'^^^ arc also not uniquely 
defined. However, this does not create any problems. 

As the interaction term Hint does not commute with 
Ho the HS transformation can be performed first sub- 
dividing the interval (0,/3) into segments of the length 
A = P/N with ^ 1 and then integrating over the 
auxiliary field for each t; = A (^ - 1/2) , ^ 1, 2...N. 

In this Section we are interested in analytical calcula- 
tions and therefore we write all formulas in the continu- 
ous limit (A — > 0) with respect to the "imaginary time" 
T. Transformations for the discrete time (finite A) will 
be performed in Section [TVl 

We rewrite the partition function Z, Eq. (|2.4|) . as 



exp ( —f3Ho ] Tr exp 



/3 



Hint (t) dr 



(2.9) 

where T^- is the time ordering operator (the time grows 
from the right to the left) and Hint (t) is the operator in 
the interaction representation. 



Hint (t) = exp (^HqtJ Hint exp (^-Horj 



(2.10) 
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The term ijf") can be decoupled as follows 



Trexp[~J^ H'^l{T)dT) (2.11 

Crcr ('^) Cr<7 

(t) dr 



M^o[0(°)] =exp(-^E/(^^°^(-))'^-) (2-12) 

where c (r) and c (r) are creation and annihilation oper- 
ators in the interaction representation [introduced anal- 
ogously to Eq. p.lOp ], (r) is a real field with the 
bosonic periodicity, (t) = (/)f°^ (t -I- /3). 

Before decoupling the term H^^l, Eq. (|2.8p . we rewrite 
it in the form 



(1) 



9 / ^ r,r' 



(1) 



2n5:0l)c+c.. + ^F(l)(2n)^ (2.13) 



where F^i) =Er' Vr'^n 



//^^!, n = (c+jCrtr) is the fermion den- 
sity per spin direction, and is the number of sites in 
the system. 

The second term in Eq. ()2.13p renormalizes once more 
the chemical potential, such that now it equals 

fi' = fi-^{Vo + V)+2nV'^'l (2.14) 

The last term in Eq. I^J^, {2nfV'^^'^Nd/2, is a trivial 
contribution to the thermodynamic potential f2. 

The fermion density n should be calculated using the 
shifted chemical potential Eq. p.l4p . and, thus solv- 
ing a self-consistency equation. Wc assume that this has 
been done and consider the first term in -fff^] in Eq. ()2.13p 
decoupling it by integration over another field (j^h^ (r). 



Tr exp 



i^i^l {r)dT 



(2.15) 



j Tr exp ( 4'' {r) {Cra (t) Cra (t) " u) dr) 



r 



where 



W, [0(1)] = exp (-\ Y: f (r) {V^'^)~l 013) (r) dr 

r,r' ^ 

(2.16) 



and 



r) is also a real periodic field, t/if^'' (r) 



</>r^'' (t + /?)■ The matrix {V^^'')^ ^, is a matrix inverse to 

the matrix V^J,. The integral over (r) in Eqs. (|2.15[ 

I2.16P converges because {V^^"^)^^, is positive definite. 
Introducing the field 

cf>rAr)^4°Hr)+'j4'Hr) (2.17) 

wc write finally the partition function Z, Eq. (|2.9p . in the 
form 



(2.18) 

with Z [(j)] equal to 



Z\ 



Tr, 



,[cxp(-/3Hc 



Tt exp 



Hl:> (r) dr 



and the Hamiltonian Hq now given by 



(2.19) 



(2.20) 



with the modified chemical potential /i' from Eq. (|2.14p . 

We calculate the trace over the fermionic operators 
c, c"*" introducing an additional variable < u < 1 , re- 
placing (pro- (t) by u(j)ra- (t) in Eq. (|2.19p and writing Z [(p] 
in the form 

Z[4>] = Zoexp[V r [ a4>ra{T) (2.21) 

L ^ ,^ Jo Jo 

X (g^!;,^} (t, t + 0) - G(°)^, (t, t + 0)) dudr 



In Eq. (j2.2ip . Zq is the partition function of the ideal 
Fermi gas. 



Zo^l[{l+cxp[-^ 



k.<7 



(2.22) 



where are the eigenvalues of the Hamiltonian Hq, 
Eq. (ESni), and 



G 



{r,r') 



TrCr. (r) (rQ exp (- /^^ H[:f dr) )^ 



TrCxp(-!^H\::^Ur)^ 

'(2.23) 

is the electron Green function in the field (t). In 
Eq. (|2.22p , the operators c, c and "f ^ (t) are written 
in the interaction representation using Hq, Eq. (|2.20p . 
as the bare Hamiltonian and G|.°^, ^ (r, r') is the Green 
function of the ideal Fermi gas described by the Hamil- 
tonian Hq. The symbol (...)q in Eq. (|2.23p stands for the 
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Gibbs averaging over the states of the Hamiltonian Hq, 

Eq. srm . 

The Green function G^""^^ (r, r') , Eq. ([2231), satisfies 
the following equation 



d_ 



[ud^ (r)] ) d^fl (r, r') = 6ry6 (r - r') , 



hrc [U(l) (r)] = Er - m' - CFu4>rcr (t) 



(2.24) 



where e^/r = ~ Tlir' ^r,r' fr' for an arbitrary function 
A conjugated equation can be written as 



d_ 



- ~ K., [u<l, (r')] ) G^^tJ. (r, r') = (r - r') , 



(2.25) 

The electron Green function G["t^ (t, t'), Eq. ([2^, sat- 
isfies the fermionic boundary conditions 

G^:,fl (r, r') = -G("^' (r + /3, r') = -G^-^) (r, r' + /3) . 



Eqs. (p:^ 1^:^ can serve as the starting 

point of our bosonization scheme. In the next subsec- 
tion, we derive equations for the bosonic excitations in 
an external time dependent HS field. 



B. Equations for bosons in the fluctuating fleld 

1. General equations. 

According to the results of the previous subsection one 
can calculate the partition function Z solving Eq. (|2.24p 
or (|2.25p for a given configuration of the field 0^0- (t), 
substitute the solution into Eq. (|2.2ip , thus obtaining the 
functional Z [0] , and then calculate Z using Eq. (|2.18p . 
Actually, solving Eqs. (|2.241 12.25^ one obtains more in- 
formation than necessary because the Green function 
G^^f).^ (t, t') that should be obtained from this equa- 
tions depends on two times, r and r', whereas these time 
should be taken equal in Eq. (|2.2ip . 

It is not difficult to derive closed equations for the 
Green function c'ff).^ (r, r -I- 0) entering Eq. ^T^. Sub- 
tracting Eq. (P?^ from Eq. ((^^ and putting t' r-t-O, 
we obtain 



d 



^ - hr„ [u<i> (r)] + hr^, [u<j, (r)] j (r, T + 0) = 

(2.26) 

We write here equal times in the functions (j)r (t) because 
we assume that they are continuous. This is not so for the 
Green functions G^^^).^ (r, r') because they have a jump 
at equal times. For discontinuous functions (f>r (t) one 
would have to take slightly different times r and t + 0, 
too. 

It follows from Eq. (P?^ that 



_9_ 



^G(';;« (r,r + 0)=0 



(2.27) 



which shows that the total number of particles in the 
system of the non-interacting fermions in the auxiliary 
HS field 4>r (t) does not depend on time. To avoid con- 
fusion, we remind the reader that we started from the 
grand canonical formulation for the model of the inter- 
acting particles, Eqs. p.l) 12. 4p . and writing about the 
particle conservation now we do not mean that we change 
the formulation of our initial model. 

For subsequent manipulations, it is convenient to in- 
troduce a new function 



Ary (z) - G^},^ (r, r + 0) - G^^";;^ (r, r + 0) (2.28) 



(«0) 



where z = (r, cr, u). 

The function Ar_r' (z) has the property 



(z) = 0, 



is periodic, 

Ary (r, cr, u) 



■ /3, u) , 



(2.29) 



(2.30) 



and. hence, describes bosons. 

With this function, the functional Z [0] takes a simple 
form 

Z [(j)] = Zq exp — / / (Jipra ij) Ar,r {z) dudr 

r.a Jo Jo 

(2.31) 

Eq. (|2.3ip is a reformulation of (|2.2ip in terms of A^^r' (z), 
Eq. (12281). 

Having written Eq. (j2.26|) we can derive a closed equa- 
tion for the function Ar.r (z) , Eq. (|2.28p . Our procedure 
is very similar to deriving the kinetic equation starting 
from equations for Green functions^. 

We rewrite Eq. ()2.26p at w = 0, and subtract it from 
Eq. (|2.26p . Using the definition of Ar^r' {z), we come to 
the following equation for this function 



d 



Mr.r' {z) = ir — Sr' ~ Ua^r,r';cr (t) 
^r.r';cr (t) = (t>ra (t) - (f>r>cr (t) 

(0) 



(2.32) 

Gj,"/, (r, T + 0) is the Fermi distribution 

tation 

(2.33) 



where riry 

function of the ideal gas in the coordinate representation 

ip(r-r') '^P 



exp {13 (e 



The solution of Eq. (|2.32p . as it stands, is not unique. 
One can easily see that there is a non-physical solution 

Ar^r' (2^) ~ Tir^r' for any (f)ra 

(t). Moreover, for static 
one can add to the solution of Eq. (|2.32[) a 



fields 



combination A^f'},, of the form 



a: 



(0) 



E 

k 



(2.34) 
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density of states at the Fermi surface 



FIG. 1: Simple electronic diagrams. 



where are eigcnfunctions of the operator hr [ucj)] and 
Cfe arc time independent coefficients. 

There can be less trivial solutions of the homogeneous 
equation corresponding to Eq. (j2.32p when (f)ra (t) is a 
discontinuous function of time with large jumps. Such 
solutions correspond to poles in the Green functions 
Gr,r' {t, t + 0) and should be treated very carefully. The 
sign problem in the MC simulations appears as a result 
of the poles in the integrals over u in Eqs. (|2.3ip which 
arise when the fields (f>ra (t) are discontinuous in time. 
We discuss this question in Section HVl 

It is clear that the ambiguity in solving Eq. (j2.32p fol- 
lowing from the existence of solutions of the homogeneous 
equation creates problems in both numerical and analyt- 
ical treatment. This means that a procedure fixing the 
proper solutions of Eq. (|2.32p is necessary and we intro- 
duce it below. 



2. Regularization of the bosonized model. 



Seeking for the unique correct solution of Eq. (|2.32p we 
are guided by the conventional diagrammatic technique 
for fermions which is well defined and all corresponding 
diagrams can be calculated at least in principle. In order 
to visualize how we come to solving Eq. (|2.32p instead 
of summing the conventional diagrams we consider the 
simplest loop represented in Fig. [1] (a). The expression 
corresponding to this diagram can be written as 



T 



ie — e'ie 

e.p I' 



luj — e: 



p+q 



E 



'p+q 



p *^-^P + q 



(2.35) 



where e' 



and 



2iT{n+l/2)T and uj 



2'KmT are fermionic and bosonic frequencies, respec- 
tively. 

The expression in the second line of Eq. (|2.35p corre- 
sponds to the solution of Eq. (j2.32p in the lowest order 
in (j), which demonstrates how the bosonic modes are ob- 
tained from the fermionic lines. 

However, the first and the second lines in Eq. (j2.35p 
are different at a; = 0, q = 0. The first line gives the 



y^cosh- (^]=?^ 

^ 2T \2T J dn 



(2.36) 



but the expression in the second line is not defined. Of 
course, one could speak about the limit q — >■ instead of 
just putting q = but this is not justified for any finite 
physical system with discrete energy levels. 

The uncertainty of the expression in the second line of 
Eq. (|2.35p is reflected in the existence of solutions [like 
those given by Eq. (|2.34p ] of the homogeneous equation 
corresponding to Eq. (|2.32p . 

We see that our bosonization scheme is not exactly 
equivalent to the initial fermionic model and needs a reg- 
ularization in order to avoid the uncertainty. This un- 
certainty appears not only in the situation represented 
in Fig. [T] (a), where the field 4> enters with w = 0, q =0. 
The same problem is encountered when calculating the 
contribution of, e.g.. Fig. [1] (b), using the bosonization 
scheme because this graph contains fermionic Green func- 
tions at coinciding momenta and frequencies but the in- 
teraction lines may carry any momenta and frequencies. 
We will come back to this point when discussing in detail 
the diagrammatics of the bosonized theory. 

The hint at a possible regularization is given by the 
form of Eq. (|2.35p . If the momenta q were continuous 
rather then discrete and we could simply neglect the con- 
tribution of the state with q =0, the bosonization scheme 
would work in this order. This is not sufhcient, though, 
because the momenta and frequencies of the two horizon- 
tal fermion lines in Fig. [1] (b) obtained after integration 
over the HS field coincide for any translationally invari- 
ant systems and one should slightly violate the transla- 
tional invariance in order to split the momenta. 

This goal can be achieved if we assume that the sys- 
tem of the interacting electrons considered here is imbed- 
ded in a "bath". We introduce this bath considering a 
model of the interacting electrons on a d-dimensional lat- 
tice with the total number of sites iV^°*'^'. However, we 
assume that the interaction Vry , Eq. p.3p , entering the 
Hamiltonian H, Eq. (|2.ip . vanishes outside a subsystem 
consisting of a considerably smaller number of sites Nd- 
As an example, we can suggest the following form of the 
interaction Vr r' 



V 



total 



Vry , |r|,|r'| <i?o 
, otherwise 



(2.37) 



which means that the interaction is finite inside the 
sphere of the radius i?o containing Nd sites and vanishes 
outside this sphere. In other words, we attach metallic 
leads to the system of the interacting electrons and as- 
sume that the inter-electron interaction vanishes in the 
leads. 

It is clear that the leads cannot change physics of the 
system of the interacting electrons if the latter is suffi- 
ciently large. However, this model is reasonable even for 
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a system with a small number sites (quantum dot), al- 
though properties of the systems with and without leads 
can be different. 

The form of the interaction Vj:p°-\ Eq. (|2?37l) . formally 
violates the translational invariance and the matrix ele- 
ments 

(q, qo) = J2 V;*°^'^'e-'i°(''+''')/2-'i("— ■') (2.38) 



are finite for qo 7^ 0. 

Of course, the matrix element with qo = remains 
finite and, for such momenta, the two Green functions 
in Fig. [T] (b) still have the same momenta and frequen- 
cies. However, we can neglect them because they give a 
small contribution provided the number of the contribut- 
ing momenta qo is large. The latter is achieved for a large 
number of the sites in the bath. 

We emphasize, that neglecting the contribution of 
graphs containing several electron Green functions with 
coinciding momenta and frequencies is possible because 
their contribution is not singular. One can estimate the 
relative error of this approximation as being proportional 
to Nd/Nf^"'-, where iV*"*"' is the total number of sites in 
the entire system including the bath. 

Following this idea we could simply exclude in the 
bosonization approach all propagators with w = 0, q = 
"by hand" . However, this is possible only when doing a 
perturbation theory and such an approach is not suffi- 
cient for non-perturbative and numerical studies. 

The bosonic modes with w = 0, q = correspond to 
solutions of the homogeneous equation in Eq. (|2.32p . In 
order to discard these modes by a regular procedure we 
regularize Eq. (|2.32|) slightly changing it. This is not a 
trivial task because this slight modification should guar- 
antee the complete absence of the solutions of the homo- 
geneous equation corresponding to Eq. p.32p . 

We regularize the bosonized model replacing Eq. p.32p 
by a system of two equations 

T-Lry [t) Aa-r.r' {z) = -uaUr.r' ,a^r.r' -^cr {t) Ba, (2.39) 

where a — 1,2, Hr^r' is a 2 x 2 matrix 

d 

Hr.r' (r) = klMr,r' {z) + iAjT^ + A7, (2.40) 

OT 



and 



Al-r,r' {z) 



M-r.r' {z) 



{z) 

^l:r,r' (^) 



Bi = 



The functions Aa-j,,r' (z) satisfy the bosonic boundary 
conditions 

Aa-r,r' (t, CT, u) = Aa:r,r' {t + P , CT, u) (2.41) 



and 2x2 Pauli matrices are used 



1 y ' ^ \i ) ' \^ -'^ 

(2.42) 

The real parameter 7 should be put to zero, 7 — >■ 0, at 
the end of calculations. 

With this modification we write the function Z [^] in 
the form 



Z [cj)] = Zq exp 



-IE 



T4>ra (r) A'^.^ ^ (z) dudr 



(2.43) 

The exponent in Eq. (j2.43p is an even function of 7. 
It is not difficult to see that putting 7 = in Eqs. (|2.39I 
I2.43P one returns to Eqs. (|2.31[ [2732)) . However, keeping 
in Eqs. (j2. 39112. the parameter 7 finite allows us to 
discard all the solutions of the homogeneous equation in 
Eq. (12321) . 

This can be understood rather easily because the op- 
erator H is hermitian and, therefore, its eigenvalues are 
real. Moreover, as will be shown in the next subsection, 
they appear in pairs: if an eigenvalue E exists, then the 
eigenvalue —E exists, too. If, when changing parameters 
of the operator T-L, an eigenvalue E with its counterpart 
—E turned to zero at some point, this would mean that 
after crossing this point they become imaginary. How- 
ever, the latter is forbidden by the hermiticity of the 
operator H. A general proof of the absence of the solu- 
tions of the homogeneous equation in Eq. (|2.39p or, in 
other words, absence of zero eigenvalues E of the opera- 
tor Tir^r' (t), Eq. ()2.40|) . is given in Appendix [A] 

The absence of zero eigenvalues makes the operator 
% invertiblc. Since all the parameters of Eqs. (|2.39I 
I2.40p are real, the solutions Aa-^.r' (z) are real, too. This 
means that the exponent in Eq. p.43p is real and Z [(p] 
is real and positive. This property of the Z [0] can be 
very important for numerical computations using the MC 
method. 

Using the regularization with the parameter 7, 
Eqs. (|2.39ll2.43p . we can prove a stronger than Eq. (|2729)) 
relation for Ar_r' (z), 



^Ar,r (z) = 



(2.44) 



Eq. (|2.44p is fulfilled automatically for any solution of 
Eqs. (|2. 39112. 40|) . This property can be proven putting in 
Eqs. (|2.39[ I2.40p r = r' and summing over r. Then, one 
obtains the equations 



^i'Az)+di':{z)/dT = 
di',{z)/dT+^i:[z) = 



(2.45) 



where (2) = ^a;.,. {z), I'^ [z) = {z). 

Only the trivial solution I'^ (z) = /" (z) — satifics the 
boundary condition, Eq. (|2.4ip . leading to Eq. (|2.44p . 
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3. Spectral expansion. 

Ill this subsection we discuss important properties of 
the solutions of Eqs. (|2.39[ I2.40p using spectral expan- 
sions in eigenfunctions of this equation. In order to better 
understand the difference between solutions of Eq. (j2.32|) 
and Eqs. (|2.391 I^IW)) . let us consider both the equations. 

The operator in the L.H.S. of Eq. ()2.32p is not her- 
mitian, which results in a double set of eigenfunctions 
{v^^, (r)} and {v^^, (r)} satisfying the equations 



^ + Mr.r' ) v^y (r) = A^<,, (r) , (2.46) 



dr 



with the boundary conditions v^^^, (r) — v^^^, (r + /?) , 

v^r' i^) — ^rr' {'^ + P) ■ Two othcr cquatious can be writ- 
ten taking the complex conjugate of Eqs. (|2.46p . 



d_ 

dr 



+ Mr,r' <; (r) - A^*<; (r) , (2.47) 



The orthogonality conditions for the functions v^^, can 
be written as 

E f<r' (r)<; (r)dr = <5^-^' 

r,r' 

(the same for the complex conjugates). 

The eigenvalues are generally complex and A^^ = 
is not excluded. For example, the function A^^^,, 
Eq. p.34p . corresponds to the zero eigenvalue for static 

All non-zero eigenvalues A^' appear in pairs. If an 
eigenvalue A^ with an eigenfunction v^^., (r) exists, then 
the eigenvalue — A^ with the eigenfunction ^ (t) ex- 
ists, too. We see that nothing prevents the eigenvalues 
from turning at some points to zero and being com- 
plex, which is an unpleasant feature of the theory. 

In contrast, the operator Jir^r' (t) is hermitian. The 
structure of the two-component vectors S^^, (r) can be 
represented as 



and they can be found from the equation 

-Hry (r) S^y (r) = E'' S^y (r) (2.48) 

supplemented by the boundary condition S^^t (r) — 
S^^, {t + (3), where the eigenvalues must be real. 

The eigenstates with eigenfunctions S^^,,, (r) and eigen- 
values E^ have their counterparts with the eigenfunc- 
tions AiS^ ^ (r) and eigenvalues —E^. 



The orthogonality condition for the eigenvectors 
S^^, (r) follows from Eq. (|2.48p and can be written as 



E 



(^^f,, (r) yS^y (r) dT = ,5^'^' (2.49) 



where the symbol "-t-" stands for the complex conjugation 
and transposition "T" . 

The completeness of the system of the eigenvectors S 
expressed by the relation 



K 



fy (r) {S^U (^i) ) ^ = ^r,rA',r',S (t - t') (2.50) 



allows us to write the spectral expansion for the 2x2 ma- 
trix Green function Gr,r';r2_,r{ {Tj'''!) satisfying the equa- 
tion 

'Hr,r' (t) Gr,r'iri,r'^ (t, Ti) = Sr^nS^y^S (t - Ti) (2.51) 

in the form 

(r) fe; (ri))"" 

gry;r,y (t, ^ ) = E " ^ (^.52) 

K 



Since all eigenvalues E^ are not equal to zero, the Green 
function Gr,r':ri,r[ {'''jT'i): Eq. (|2.52p . is uniquely defined 
and does not have singularities. Using the properties of 
the eigenvectors S^^, (t) explained after Eq. ()2.48p we 
obtain a symmetry property for these functions 

Qr,r';ri,r[ (TjTi) = - A-lQr' ,r;r[,ri (tIjT) Ai (2.53) 



With the Green function Gr.r':ri,r[ {t,ti) , Eqs. (|2.51[ 
I2.52p . we write the function Z [0], Eq. (j2.43p . in the form 



Z [(j)] = Zq exp 



E 



(l)ra (r) 



(r) drdu 



(2.54) 



Eqs. (j2.52l I2.54P give an explicit unambiguous expres- 
sion for the function Z [(j)] . It can be further integrated 
over 0rcr using superfields, as it is done in the next Sec- 
tion, or calculated numerically. In the latter case, it is 
more convenient to solve directly Eq. (|2.39p and substi- 
tute the solution into Eq. (|2.43p . 

We check our regularization scheme in Appendix IB] on 
the model with a static HS field (f>rcr- 

One can interpret the real function Ar^r' (z) as a den- 
sity matrix, while Eqs. (|2.32l 12.39^ are analogues of the 
von Neumann equation. In this language, Eq. (|2.44p is 
the particle conservation law. Using the density matrix 
in quantum mechanics implies a connection of the system 
to an environment. Within our approach, we use the reg- 
ularization based on the introduction of the bath. This 
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enforces in a very natural way the analogy between the 
matrix A and the density matrix in quantum mechanics. 
One can say that our bosonization procedure corresponds 
to replacing the description of the quantum mechanics in 
terms of the wave functions by the description in terms 
of the density matrix. 

In the classical limit, Eqs. (j2.32[[2?39|) play a role of the 
kinetic equation and (z) (after Fourier transforming 
in r — r') can be considered as the distribution function 
for a particle moving in a fluctuating field in imaginary 
time. 

We see that the original problem of calculating the 
partition function of interacting fermions Z, Eq. (j2.4p . 
has been reduced after the rcgularization to solving 
Eqs. (|2.391 I2.40p for a bosonic function Ar^r' {z), sub- 
stituting the solution into Eq. (|2.43p and calculating the 
functional integral over (j>rcr (t) in Eq. (|2.18p . For an an- 
alytical investigation, this integral over (j)ra- (t) can be 
performed before doing any approximations and this is 
demonstrated in the next Section. 



III. SUPERFIELD THEORY OF INTERACTING 
BOSONS AND PERTURBATION THEORY 

A. Superfleld theory of interacting bosons 

In this subsection, we reduce the computation of the 
partition function Z, Eq. (|2.4p . to calculations with a su- 
perfleld model describing interacting bosonic excitations. 
Eq. (12.321) resembles quasiclassical equations written in 
Ref. [la- The solution of these equations was represented 
in terms of a functional integral over 48-component su- 
pervectors, which allowed the authors to integrate over 
the HS field. Unfortunately, the resulting Lagrangian 
and calculations with it were rather cumbersome due to 
the large number of the components of the supervectors. 

Now we use another trick, known as the Becchi-Rouet- 
Stora-Tuytin (BRST) transformation^^ (see also the 
book by Zinn- Justin^) . A similar transformation was 
used in the quantization of non-abelian gauge theories^. 
In condensed matter physics, this trick has been used 
for the first time in Ref. [s^ Using this transformation 
one can also represent the solution of an equation (not 
necessarily linear) in terms of a functional integral over 
both conventional and Grassmann anticommuting fields. 
Due to some special symmetries, both the types of the 
fields can be unified in a superfields depending not only 
on coordinates and times but also on additional anticom- 
muting variables. 

In order to avoid complicated formulas containing 2x2 
matrices that appear as a result of the rcgularization, we 
restrict our consideration in this Section by the thermo- 
dynamic limit. In this limit the spectrum is continuous 
and the states with coinciding momenta and frequen- 
cies do not give any essential contribution anyway. This 
means that we put the parameter 7 = 0. At the same 
time, we still assume that the system is imbedded in the 



bath. This allows us to treat irregular expressions in a 
simple way because the translational invariance of the 
subsystem (where the interaction is finite) is broken and 
Green functions with coinciding momenta and frequen- 
cies give a vanishing contribution in the thermodynamic 
limit considered in this Section. 

Let us explain this procedure in more details. Suppose, 
we have an equation 



F(A)-0, 



(3.1) 



where F (A) is a real matrix function of a real matrix 
function A, and our task is to find the solution Aq of this 
equation and calculate a quantity B (Aq) , where B (A) 
is another matrix function. Instead of proceeding in this 
way, one can write B (Aq) in a form of the integral 



BiAo) = J B{a)S[F{a)] 



det ( — 

oa 



da (3.2) 



over the variable a of the same structure as A. The mod- 
ulus of the determinant in Eq. (j3.2[) is included in the 
integrand in order to normalize the (5-function. 
As the next step, we represent the (5-function as 

6 [F {a)] = C / cxp [iTr {fF (a))] df (3.3) 



where / is a real matrix variable and the symbol "T" 
means the transposition of r and r', /J^, = fr'.r . C is a 
normalization constant. 

The determinant in Eq. p.2|) can also be represented in 
a form of an integral but now the variables of integration 
should be anticommuting Grassmann variables rj and 77+ , 



det 



dF_ 

da 



exp 



1 Q^^ 



d7]+dri. (3.4) 



After unifying the fields a, 77, and 77"*" in a super- 
field 4*, the expression in the exponential can be written 
in a compact way. We do not present here general for- 
mulas that can be found in the book [3l| but concentrate 
exphcitly on Eqs. (|2.31l2.32p . 

Eq. (|2.32p is linear and real and we can apply Eq. (|3.2p 
directly. Using the transformation of Eq. p.3p we obtain 
in the exponent a quadratic form in terms of the variables 
a and /. Then, we write the determinant with the help 
of Eq. p.4p and obtain another quadratic form in the 
exponent containing the anticommuting variables rj and 

At first glance, there should be a problem related to 
the presence of the modulus in Eq. p.2p and absence of it 
in Eq. p.4p . Fortunately, the operator d/dr + Mr.r' (z) is 
real and antisymmetric, which leads to an always positive 
determinant. Therefore, the modulus in Eq. p.2p does 
not play any role in the case under consideration. 

Writing the quadratic forms in the exponent using the 
variables a, /"^, 77, and 77"*" is not difficult and we do 
not write them in this form. A more compact form of 
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the functional integrals can be achieved introducing the 
supcrficlds (R), R = {z,d,9*}, as follows 

^r.r' (R) - ar.r' {z) 9 + fj^, (z) 9* + Tj^y (z) + 7?+^, {z) 9* 9 

(3.5) 

where a^y {z) and f^.^ (z) are real fields of the coordi- 
nates r, r' and the variables z = (r, cr, u) . The hermi- 
tian conjugation means complex conjugation " * " sup- 
plemented by the transposition, such that r/^^, = vt.'r- 
The variables 9, 9* are artificially introduced Grassmann 
variables that help us to write the exponent in the com- 
pact form. 

As a result, we write the functional Z [(f)], Eq. (|2.3ip . 
in the form 

Z [(f)] = Zo J exp (-5o - [*]) (3.6) 



where So is the bare part of the action 

d 



(3.7) 

The second term ^S^""^) [^f] in the exponent in Eq. ([3J)) 
is linear with respect to the HS field 0rcr (t) ■ Its explicit 
form reads 

[^] = -iY,j (t>r. (r) [7/ (*.',. (i?) - rv'^r^) 

r,r' 

X (*^^r, (i?) - nr,r'0) - i(T<5^,^/5'r,^ (i?) 61*1 dR (3.8) 



It is interesting to remark here that the superfield ^ (i?) 
is anticommuting, which is a rather unusual feature 
of the field theory considered here. However, it is 
very important that this field describes bosons and not 
fermions. This follows from the periodic boundary con- 
dition (r) = ^E" (r -I- /3) . It is this condition that deter- 
mines unambiguously the statistics of the particles. The 
fact that (i?) is anticommuting is only a formal prop- 
erty 

The form of Eqs. (|3.6ll3.8p allows us to average im- 
mediately over the HS field (pra {t). This integration is 
Gaussian and is specified by Eq. (|2.18p . The analytic 
supersymmetric field theory written here relies crucially 
on the presence of a bath, whose function is to break the 
translational symmetry of the system. In order to take 
the bath into account, we model the interaction in the 
same way as in Eq. (|2.37p . The pair correlation of the 
fields with the distribution W can be written as 

{aa'cj^ra (r) c^',,„, {t'))^- = Ury (i?, R') (3.9) 

Ur,r' {R, R!) =5{t- t') (aa'V^°''''5r,r' + K^'''*"*"') - 



where Vq*"*"' and ^ 



{l),total 



are defined from Vq and y*-^^ 
according to Eqn. (|2.37p . The coupling constants Vo and 
y(^) contain the constant V that has been introduced 



in a rather arbitrary manner (see, below Eq. (|2.5p . How- 
ever, the same constant enters the renormalized chemical 
potential Eq. (|2.14p . and the renormahzed thermody- 
namical potential Vl written below Eq. (|2.14p . Of course, 
the constant V should disappear from the final result 
for the thcrmodynamical potential, which can serve as 
a check of any computation. Actually, considering only 
perturbation theory in the interaction there is no neces- 
sity to introduce this constant. In such calculations the 
sign of Vry in Eq. (|2.4p does not play an important role 
because one could decouple the interaction term with V 
integrating over purely imaginary HS. We have added V 
because we want to keep the HS fields real, keeping in 
mind a possibility of applying the method to numerical 
investigations. 

Using Eqs. p.8l3.9p we easily integrate in Eq. (|3.6p 
over the field (p^u (r) reducing the partition function Z 
to the form 

Z = ZoJ exp {-So m - S„,t [*]) (3.10) 

with So [*] given by Eq. ([XT)) and 

5™t [*] ^ S2 [*] + 53 [*] + 54 [*] . (3.11) 

The term Sint [^] h^ the action describes the interac- 
tion between the fields The terms ^2 [^], S3 [^] and 
54 [^'] contain quadratic, cubic and quartic in 4* terms, 
respectively They can be written in the form, 

S2 = \jY. {^'^tARW - {R) , fl]rrU9} Ur^r, [R, Rl) 
r,ri 

X {l^r,,rARlWl - {Rl) .n\r^,r,Ui9^} dRdRi, 

(3.12a) 

S3= J *r',.(i?)*r,r'(i?)f/.,ri(i?,fll) 

r.r' ,ri 

X {i'^ri.ri{Rl)Ol - {Rl) ,n]r^,r^Ui9i} udRdRi, 

(3.12b) 

Si^\j J2 '^r'AR)'^r,r'{R)Ur.rAR,Rl) 



r,r ,ri .r-^/ 

X ^r,^ ,.^{Ri)^r^ ,.,^{Ri)uiudRdRi 



(3.12c) 



where [, ] stands for the commutator. Integration over R 
in Eq. p.l2all3.12cp implies summation over a and inte- 
gration over u, T, 9, 9* . The bare action So and the inter- 
action term 54, being invariant under the transformation 
of the fields ^ 



+ K*) 



(3.13) 



{k and K* being anticommuting variables) are fully su- 
persymmetric in the sense of Ref. Isil . whereas the terms 
52 and ^3 break this invariance. The invariance under 
the transformation ()3.13p is stronger than the standard 
BRST symmetry for stochastic field equations [invariance 
under the transformation {9*) '9 {9* + k*)]^, and 
reflects additional symmetries of Eqs. ()3.7l3.12cp . 
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Eqs. (| 3. ion 3. 12c)) completely determine the new 
bosonic superfield theory. This model can be studied 
using standard methods of field theory. One can, e.g., 
expand in the interaction U using the Wick theorem with 
simple contraction rules following from the form of the 
bare action Sq. In the next subsection, we will demon- 
strate how such calculations can be carried out explicitly 
but now let us understand what one obtains neglecting 
the cubic and quartic in ^ terms ^3 and ^4 in the action. 
In this approximation, one has a purely quadratic action 
and, making the Fourier transform with respect to time 
and space, we obtain easily for the partition function Z 
an RPA-like expression. 



Z ~ Zq exp 



2 (27r) 

'^p-k/2 - «p+k/2 



(3.14) 



ILjJ - 



-p-k/2 



-k/2 {2TTf 



The same result can be obtained using Eqs. (|2.18|2.3H 
I2.32[) and neglecting the field 4>r(j{T) in the L.H.S. of 
Eq. (12321) . 

In Eq. p.l4|) . {K — 1) is the contribution of non- 
interacting bosonic excitations. Considering their inter- 
action, one can fully describe the initial fermionic system. 

We note that the boson model, as formulated in this 
paper, does not appear convenient for a standard analyt- 
ical approach such as renormalization group procedures, 
where the low lying modes should have been singled out 
in the very beginning. However, the superfield model 
Eq. p.ll|) and a proper low-energy formulation arc for- 
mally quite similar. Therefore, this Section does not only 
serve as check of the superfield model but also as a guide 
for future analytical studies of low energy effects. 

In the next subsections, we demonstrate how the con- 
tributions of the conventional perturbation theory^ in the 
original fermion language are reproduced in the boson 
superfield representation. We perform this check calcu- 
lating the perturbation series for the thermodynamic po- 
tential Q up to the second order in the interaction. 



B. Wick theorem and diagrams 

In order to develop the perturbation theory in the in- 
teraction, we Fourier transform the superfields 

where p = (cr, w, 6, 9*) and w stands for the bosonic Mat- 
subara frequencies. The momentum integration (dp) = 
dp/(27r)'' extends over the first Brillouin zone. 

In the Fourier representation, the bare action reads 

'^0 = / *P'P.--(P){*^ - [ep - £p']}'fpp'Ap) 

X {dpdp')dp (3.15) 



and similarly one can obtain the Fourier transform of the 
interaction, Eqs. p.l2all3.12cl) . 

According to Eq. p.isp . the bare propagator of the 
bosonized theory is given by 

{^p,p'■Ap)^PuP[■,^^ (pi)) ^{0- OiW - ei)5{u - u,) 

X 5.., {.2nf''S{p - p[)6{p' - p,)^-^g{uj- p, p') 

(3.16) 

wherein, g{uj; p, p') is the Green function in the momen- 
tum representation 



ff(w;p,p') = i [i 



Using Wick's theorem and the pairing rule, Eq. p.l6p . 
one can calculate averages of \l/-products arising in per- 
turbation theory expansions analytically in a standard 
manner. 

For a convenient diagrammatical representation, the 
propagator can be depicted in the form of an oriented 
double-line as shown in Fig. [2] (a), the first momentum 
p being carried by the right-facing upper line, the second 
momentum p' by the left-facing lower line. 

In the same spirit, the interaction vertices, Eqs. (j3.12al 
I3.12c[) . are diagrammatically represented as shown in 
Fig. O (b). Depending on whether a leg is depicted as 
"out-going" [right-facing] or "incoming" [left- facing] , the 
momenta and frequencies are either reversed [momenta p 
and p' arc interchanged and frequency co changes the 
sign] or not. The interaction lines carrying the momen- 
tum q correspond to the propagator 



yd) 



(3.17) 



and are attached according to the analytical structure of 
Eqs. p.l2all3.12cp . In this way, the single lines within 
the double-line structure obey evident momentum con- 
servation rules, whereas the bosonic frequency of each 
double-line is conserved in the usual way and, therefore, 
not labelled in the diagram. 

With the rules formulated, we are now ready to discuss 
the perturbation theory for the thermodynamic poten- 
tial ri. 



C. First order of the perturbation theory 

As 54 does not break the symmetry specified by 
Eq. p.l3p . the first order contribution to il in the su- 
perfield theory comes only from ^2. However, T{S2) 
is not the only contribution to the first order correc- 
tion Arid-* in the interaction potential. In addition, 
one should renormalize the chemical potential thus us- 
ing /i', Eq. (|2.14p . and shift the thermodynamic po- 
tential 51 ri -I- (2n) yd-'A"d/2 in accordance with 
Eq. (|2.13p . This leads to a trivial first order contribution 
to the thermodynamic potential 51 that can be written 
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(b) 




p + q 



p + q 



p + q p' p' p' p' 

FIG. 2: Diagrammatic building blocks for perturbative cal- 
culations: (a) the boson propagator g{u};p,p'), Eq. (|3.16|l . 
(b) the interaction terms, Eqs. (|3.12al3.12c|l . 




FIG. 3: First order diagram from (52 



as 



(5f7(i) = Nd[2n'^V'^^^ (3.18) 

where V'-^^ = - J2r' ^ry + V and n is the (unperturbed) 
fermion density per one spin direction at site r. The first 
term in Eq. (j3.18p comes from the shift of the thermody- 
namic potential, Eq. (|2.13|) . The second line originates 
from the shift of the chemical potential, Eq. (|2.14p . and 
is obtained using the standard relation dfl/dii — ^2nNd 
(factor 2 is due to spin). So, the total first order contri- 
bution Afi^^^ to the thermodynamic potential takes the 
form 



Ai7(i) = T{S2) + sn^'^ 



(3.19) 



Diagrammatically, T{S2) is represented in Fig. |31 The 
correspondence of Fig. [3] to the conventional diagram- 
matic technique becomes evident if we associate the 
double-line with two single lines representing the origi- 
nal fermionic propagators. Then, one can easily imagine 
that the diagram in Fig. [3] should yield the Fock contri- 
bution [Fig. [T] (a)]. Actually, this is not completely cor- 
rect because this diagram contains also some other terms. 
Calculating T{S2) with the help of the contraction rule 



(|3.16p and Eq. p.l7|) we obtain in the first order 

T{S2)^T^Y. I (*P,P+q;-(p)*P'+q,P';-(Pi)) 
X {-i)uiUq.^{a,(Ti)9*9i[nl,^^- n^,] dpdpi{dpdp'dq) 



-p+q 



(3.20) 



where rip is the bare Fermi distribution. 

We can represent Eq. p.20p in a somewhat different 
form using the identity 



"p' "p 



1 



■,0 



coth ■ 



p' 



2T 
(3.21) 



Substituting Eq. ([33T|) into Eq. ((3?20| . we see that the 
term with Vq in Eq. p.20p , indeed, participates in form- 
ing the Fock contribution Afi^'' to the thermodynamic 
potential 

AfiW = y VqnOnO+q(dpdq) (3.22) 

However, there are other contributions STl'-^'' in 
Eq. ((3?20l) that should be added to that of Eq. (|3T8)) 

^NdJ [2n;^nO+q (K,, + V) 
+ {nl + nO^q) (yq/2 - Vr,r ~ V) ] (dpdq)(3.23) 



The structure of the integrand in Eq. p.23p allows one 
easily integrate over p and q. Then, adding Eqs. ()3.181 
I3.23P to each other and using Eqs. (|2.5ll2.7p we obtain 

AQ^jP ^ jrj(i' + = 2Ndn^Vq=o (3.24) 

which is the standard Hartree contribution. The artificial 
interaction V introduced in the model drops out from the 
final formulas as it should. 

Thus, we come in the first order to the standard ex- 
pression for the correction Afi'^^ to the thermodynamic 
potential 



(3.25) 



where A17^^ and Afl'-p^ are the Hartree, Eq. ([3^ . and 
Fock, Eq. (|3.22p . contributions, respectively. 



D. Second order of the perturbation theory 

The perturbation expansion in the conventional field 
theory for interacting fermions leads to the diagrams 
shown in Fig. S) Only the diagrams in Fig. U] (a) and (b) 
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FIG. 4: Diagrams for the second order in the fermion language. Interesting low-energy contributions are only due to diagrams (a) 
and (b) whereas diagram (c) just renormalize the chemical potential. To check the equivalence of the bosonized model to the 
original fermion one, we must check, however, that all fermion diagrammatic contributions are exactly reproduced. 



are important for low-energy physics, whereas the dia- 
gram (c) merely modifies the chemical potential. The 
diagrams 2] (d.l) and (d.2) contain the Hartree bubbles 
that also result in a renormalization of the chemical po- 
tential. 

Nevertheless, in order to understand details it is im- 
portant to obtain all contributions of the perturbation 
theory for the boson superfield theory. Checking the ex- 
act correspondence between the original fermion and the 
boson models is what this Section is devoted to. 

Before presenting the details of the calculations, let us 
make some general remarks. For ijf^] [Eq. ()2.8p ]. the 
Hartree-type contributions H] (d.l) and S] (d.2) are in the 
boson theory generally accounted for by the renormal- 
ization of the chemical potential [Eq. ()2.14p ] and by the 
trivial contribution to the thermodynamic potential in 
Eq. (|2.13p . There is no contribution of diagrams con- 
taining closed loops of boson propagators because they 
vanish due to the superfield symmetry. The effective sec- 
ond order diagrams in the bosonized theory arc shown in 

Fig.m 

Drawing parallels with the conventional fermion dia- 
grams we interpret the double-lines in Fig. [5] as pairs of 
single fermion propagator lines. This is similar to what 
we successfully did when discussing the first order dia- 
grams. In this spirit, inspection of the diagrams sug- 
gests that the conventional RPA diagram Fig. |4] (a) is 
reproduced in the boson language by Fig. [5] (a) , the cor- 
relation diagram Fig. H] (b) by Figs. [S] (b.4,c.l), and, fi- 
nally, the conventional diagram with the two Fock loops. 
Fig. m (c), by Figs. [5] (b.l-b.3,c.2). However, for con- 
tributions coming from the on-site interaction Hamilto- 



nian H 



(0) 



Eq. 



this "graphical picture" describes 



the correspondence in a less strict way since the conven- 
tional diagrams (a) and (b) as well as the conventional 
diagrams (c), (d.l), and (d.2) in Fig. |3] coincide in this 
case, thus making a parallel between just Fig. H] (a) and 
Fig. [5] (a) meaningless. In this case one can speak only 
about a "weaker" correspondence between Figs. U] (a,b) 
and Figs. [S](a,b.4,c.l,c.2). 

Having established these parallels, the detailed calcu- 
lations can be performed in an organized manner, which 
is straightforward though a little tedious. 

Starting with (^2 ), the evaluation of the diagram in 



Figs. E] (a) yields 

X [np'+q - np^]g{uj, p + q, p)5(w, p' + q, p'){dpdp' dq), 

ee'uj 

X GO ,p,+qGO+„,p,(dpdp'dq) (3.26a) 



xG°+ 



/(dpdqdq') 



(3.26b) 



where Eq. ((235|) 

^ i^e - Cp Y 



^e.p 



has been applied and the notation 
^ for the fermion Green fmiction in 



the original electron language is used, $p = £p — 
the second term (|3.26bp , the momenta in the Green func- 
tions have been shifted using the fact that at least one 
of the involved interaction couplings V is independent of 
the running momentum q (on-site interaction). 

Inspection of the expressions obtained shows that the 
first term (|3.26al) is indeed identical to the conventional 
second order RPA contribution. Fig. |3](a). This one-to- 
one correspondence of the Vq-term confirms the valid- 
ity of the "graphical correspondence" . The form of the 
second term resembles the contribution of the diagram 
Fig. m (b) and therefore we expect that it should be can- 
celled by those diagrams of {83'^) and (52^4) that we 
have already related to the conventional contribution of 
Fig. 21 (b). Since for the purely on-site interaction the 
diagrams Fig. |4] (a) and (b) are indistinguishable, there 
is no surprise that (^2^) also contributes to the conven- 
tional diagram in Fig. |3] (b) due to the presence of the 
-channel, Eq. (1^ . 

So, let us consider now the boson diagrams shown 
in Figs. [5] (b.4,c.l) which are of conventional type of 
Fig. m (b) in terms of their graphical shape. Their eval- 
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(b.l) 



(a) 



S2 S2 




(b.3) 



FIG. 5: Second order diagrams describing perturbation series of the superfield theory. Contributions come from {S2^) [(a)], 
{S:^) [(b.l)-(b.4)], and (52^4) [(c.l),(c.2)]. 



uation yields 



{Si 



>(b.4,c.l) ^ 



X iup g{LO + io' ,p + q + q' ,p) 

X 9{oJ, P + q', p),9(w', p + q, p) (dpdqdq')- (3.27) 



Recasting Eq. p.27p into the original fermion language, 
we can equivalently write 



(Si 



,(b.4,c.l) 



-N, 



idq' 

(3.28) 



X ^£,p'-^e+w,p+q^£+aj',p+q''-^e+w+aj',p+q+q' (c'po'qC^q )• 



Inspecting Eq. (|3.28p . we see that the term with 
4:Vq — 4:VoVq exactly cancels the contribution (j3.26bp . 
while the term with V^jV^j' in Eq. (j3.28p is exactly the con- 
tribution originating from the conventional fcrmionic di- 
agram in Fig.|4](b). Thus, we conclude that we have ver- 
ified the one-to-one correspondence of the conventional 
diagram family given by Fig. 3] (a,b) and their graphical 
counterparts in the boson diagrammatics. 

Now we address the remaining second order diagrams. 
The diagrams of Figs. [5] (b.l) and (b.2) involve a formal 
propagator g(0; p, p) ~ 1/0 and we need a regularization 
scheme. This badly defined boson propagator appears as 
a counterpart to the conventional diagrammatic struc- 
ture of two fermion lines which carry identical frequency 
and momentum, cf. the discussion in Section [ll Bl How- 
ever, as has been discussed previously, this irregular ex- 
pression is healed as soon as the interaction Hamiltonian 
is no longer translationally invariant. 

As the main purpose of this Section is to check the per- 
turbation series, it is sufficient to perform the regulariza- 



tion by smearing the (5-function in the interaction matrix 
element V^Sqq'- By this procedure, the momenta and fre- 
quencies in the conventional diagram Fig. |4] (b) become 
in general mutually different, while the central propa- 
gator in the boson diagrams in Figs. [5] (b.l) and (b.2) 
has two generally different momentum arguments, e.g. 
g(0;p,p + q' — q). The smearing of the (5-function is 
achieved by breaking the translational symmetry of the 
interaction introducing the bath as it has been done in 
the preceding Section. 

Writing the proper analytical expressions wc obtain a 
regularized boson propagator (/(O; p, p -I- q' — q) ~ 1/ Ae, 
where £p_|_q'_q ~ £p-|-Ae with small Ae. Then, wc should 
expand in Ae all other terms. As the original fermionic 
contributions are not singular in the limit Ae — > 0, we 
expect that the small parameter Ae should be compen- 
sated. 

Calculating the diagrams in Fig.[5](b.l) and (b.2) with 
the help of this regularization scheme, then evaluating 
the remaining diagrams (b.3) and (c.2) that are regular 
by themselves and, finally, combining all these contribu- 
tions we come to the expression 



-2iVrfT2 J (214) - Vq)i2Vo - Vq,) 



X {« «p+q a'^i^; p, p + q).9(^^'; p + q', p + q) 

+j np+q. g^{uj'; p, p + q')g{uj; p + q, p + q') 
-i ?ip g'^iuj; p + q, p)g{u}': p + q', p) 
~i rip 5(0;; p + q, p).9^(w'; p + q', p) 

g{uj;p + q,p)g{uj';p + q' ,p)j (dpdqdq'). 

(3.29) 



5/x 



We sum over the Matsubara frequencies and obtain using 



15 



transformations similar to those leading to Eq. (|3.2H) 

^52^^^(b.i-3,c.2) = I (2^0 - V^){2V, - 

X np(np - l)(2np+q - l)(2np+q/ - l)(dpfiqdq') 

(dpdqdq')|V^ql/q'^np+q?lp+q, (3.30a) 

+ {W^ - 2VoVc, - 2FoVq0^np+qnp+q, (3.30b) 

- -(2Vb - V^){2Vo - ■(^q')^(2rip+q + 2?lp + q' - 1). 

(3.30c) 

It is clear that all singularities have disappeared. An- 
alyzing the expressions above we find that the first sum- 
mand (|3.30ap equals the contribution of the conventional 
diagram Fig. |4] (c). Thus, we have so far succeeded to 
show how the boson model reproduces all the conven- 
tional diagrams without the Hartree bubbles. 

Finally, we collect all trivial second order contribu- 
tions arising from the shift of the chemical potential, 
Eq. p.l4p . and the thermodynamic potential, Eq. ()2.13p . 
Note that in this procedure, we have to expand also the 
Fermi functions in the contribution of the first order di- 
agram, Eq. (|3.20p . Combining these trivial second order 
contributions with the contributions (j3.30bp and (|3.30cp 
yields the analytical expression for the collection of the 
Hartree-type diagrams Fig. U (d.l) and (d.2). The cal- 
culation is completely analogous to the one for the first 
order diagrams but is a little bit more tedious. 

Thus, the calculations performed in this Section helped 
us to see that the boson model yields up to the second or- 
der the same perturbation series as the conventional the- 
ory. The perturbative calculations in the boson model are 
somewhat more cumbersome than the standard pertur- 
bation theory. However, our motivation for constructing 
the boson model is to develop a new tool for studying 
both analytically and numerically the low temperature 
physics, rather than considering lowest order diagrams 
exactly. The calculations presented in this Section are 
merely a check of the bosonic approach. 

We note that the rcgularization based on the break- 
ing of translational symmetry is crucial only for a class 
of diagrams which usually do not play an important role 
in the framework of analytical investigation of low lying 
excitations. Therefore, studying, e.g., the non-analytical 
corrections to the specific heat of an interacting Fermi 
gasi^iSij the inherent ambiguity of Eq. (|2.32p is rather 
harmless because the important contributions do not re- 
quire any additional rcgularization. 



IV. DISCRETE TIME REPRESENTATION AND 
MONTE CARLO SIMULATIONS 

A. Quantum Monte Carlo 

Using exact transformations, we have mapped in the 
previous Sections the original fermion model, Eqs. (|2.1f 
12. 4p . onto the boson model, Eqs. (|2. 39112. 4^ . or Eqs. 
p. 10113.1^ . We have shown how one can treat the latter 
model diagrammatically making expansion in the inter- 
action terms and demonstrated in the first two orders in 
the interaction that the perturbation theory agrees with 
the conventional perturbation theory for the fermions. 

In this Section we want to compare the two equiva- 
lent models from the point of view of their suitability 
for Monte Carlo simulations. It is well known22ri2^ that 
studying models for fermions one encounters the nega- 
tive sign problem and the main advantages of the MC get 
lost. At the same time, the sign in the MC sampling is 
always positive for bosonic models. Therefore, it is quite 
natural to investigate the suitability of the boson model, 
Eqs. (|2.39II2.43P or (|3. 10113. 12cp . to MC simulations and 
clarify the question about the sign problem. 

Before starting our investigation of the problem we 
would like to emphasize that we do not attempt dis- 
cussing here effectiveness of various MC algorithms or 
compare them with each other. Our goal is to under- 
stand semi-quantitatively where the problems come from, 
how one can overcome them and how one can do com- 
putations in principle. We are aware of the fact that 
the concrete schemes of the calculations we suggest at 
the end may be not necessarily the most efficient ones, 
but writing explicit algorithms is beyond the scope of the 
present paper. 

The MC method is used for computation of the average 
of a quantity A over configurations c in a space of these 
configurations 

(A) = Z-i5^^(c)p(c), Z = Y,p{c) (4.1) 

c c 

According to this approach one chooses a set of M 
configurations {c^}, i ~ 1,2, ...A/ with the probabil- 
ity Z^^p (ci). Then, instead of calculating the entire sum 
in Eq. (|4.ip . one computes the sample average 

M 

This approach is very efficient because the time of a nu- 
merical computation grows with the size of the sample in 
a power law, which contrasts exponential growth when 
using other methods. 

Of course, such an approximation is possible only pro- 
vided all p (c) are not negative. In some cases the proper 
p (c) can be both positive and negative and formally one 
may not use the approximation (|4.2p . In this case one 
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may use instead, e.g., |p(c)| as the probability rewriting 
Eq. as 



J:,A{c)s{c)\p{c)\ fj:cHc)\p{c)\ 



where s (c) — sign[p (c)]. 

One can perform the computations starting with 
Eq. (|4.3p but the average sign s (c) becomes in many 
important cases very small, which leads to exponentially 
large computation time. 

The traditional scheme22r2^of the quantum MC simu- 
lations for the fermion models is based on the HS trans- 
formation and using Z [</>], Eq. (|2.19[) . as the probability 
p (c) for the MC sampling. The problem arising in this 
approach is that for some important configurations of 
(f'rcr (t) the functional Z [(/)] becomes negative. There arc 
some special cases like models with a half-filled band or 
with attraction between the fermions when Z [(j)] is neces- 
sarily positive. However, generically, the MC simulations 
are inefficient for the fermionic models. 

In principle, one could propose just to use the 
bosonized model, Eqs. p. 10113. 12c)) . However, it contains 
the superficlds and one cannot apply the MC method 
immediately. This diffiuculty can be overcome by decou- 
pling the interaction with the help of HS transformation. 
Of course, proceeding in this way we would come back 
to Eqs. (|2.39I - [^I151) . Now we can argue that the solution 
Ar.r' (z) of Eq. (|2.39p is necessarily real and substituting 
it into Eq. (|2.43p one must obtain positive values of Z [(/)] . 
We have explained how one comes to Eqs. (|2. 39112. 43|) us- 
ing the regularization based on introducing a bath and 
the small parameter 7. All the transformations leading 
to (|2. 39112. 43]) were justified and therefore one may enjoy 
computing using the scheme free of the sign problem. 

We think that this way of reasoning is correct, but it 
is definitely somewhat artificial: one introduces super- 
fields, integrates over the HS fields but then returns to 
the equations (|2. 39112. 45)) . After making these transfor- 
mations the sign problem desappears but why this has 
happened is not clear. 

So, it is very instructive to try to understand in a more 
direct way how the sign problem is overcome, what arc 
the approximations made, etc. Of course, one should 
discuss these things without using the superficlds. 

This is what the next subsections are devoted to. 



B. Discrete time representation 

We start the discussion with a slight extension of the 
original fermionic model, Eqs. (|2.1j|2.4p . The interaction 
between the electrons Vry, Eq. (|2.3p . does not depend on 
time. As a result, fluctuations of the HS field cpra (t) are 
not correlated in time, Eq. p.9p . The function S {t — t') 
describing these correlations does not create difficulties in 
analytic calculations but it is more convenient to slightly 



smear it for our dicussion. So, we change the model re- 
placing the weight Wq {4>) in Eq. (|2.12p by 



(4.3) [0] = exp 



1 

2Vo 



E 



i2, ^ , 2 1 C><f>r (t) n - 



(4.4) 

To simplify the discussion we neglect the interaction V^'^'' 
in Eq. (|2.8p thus restricting ourselves with the on-site 
repulsion. According to the accepted regularization with 
the help of the bath we assume that the interaction is 
equal to Vq inside the sphere with the radius Rq but 
vanishes outside this region (c.f. Eq. (|2.37p '). We are 
interested in the limit a — ?■ but let us take this limit 
at the end of the calculations. We do not see any reason 
why this limit should lead to results different from those 
for the original model. 

The presence of the derivative in Eq. (|4.4p introduces 
correlations between <f>r (r) at different times. We are 
mostly interested in the low temperature limit when the 
sign problem is strongest. In this limit, the pair corre- 
lation function between the fields 0^ (t) inside the inter- 
acting region takes a simple form 

/A / \ A f T/A exp(- \T-T'\/a) 

{(pr (r) (pr' (T )) = Vodr,r' (4.5) 

Za 

showing the time correlations on the scale a. 

In order to proceed with the MC computations, one 
should discretize the time and basic equations, which 
can be done in many different ways. The standard pro- 
cedure is to subdivide the interval (0,/3) into slices of 
the length A = p/N <C 1 and consider discrete times 
Ti = A{1 — 1/2). Within this scheme one writes the par- 
tition function Z in the form 



Z = lim 

A-i-O 



fZ' [4>]W:,[c^]Dct> 

jw;,[cp]D<j> ' 



(4.6) 



where 



Z'\ 



Or, I 



N 

exp (-Hq/t) n '^^P E 

1=1 r,a 



'r,lCr,lCr,l 



H':i*i ^»p(-^i:i:(*;. 

r 1=1 

+ {Ot/^f (0r,/+l/2 - 0r,/-l/2)^)), 



and the subscript I of (j) and c means that these functions 
are taken at the times rj. 

The disadvantage of Eq. (|4.6p is that one cannot bring 
the trace of the fermionic operators c, c into an useful 
form without making additional approximations. There- 
fore, we use a slightly different representation that allows 
us to calculate this trace exactly and derive the bosonized 
form. 
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t/A 



FIG. 6: Typical form of the function (f>r (t) 



Instead of using Eq. (|4.6p we write the partition func- 
tion Z as 



Z — lini 



J Zf[4>]w[4>]D4> 



A^o JW[(f>]D<f) 



(4.7) 



where (f)^ (r) = 0^^; for — 1) A < r < and the func- 
tional Z equals 

=Tre,c[exp(-/37?o) (4.8) 

xTrBXpi- y2<^4>rcr{T)Crcr{T)Cra{T)dT] , 

The weight W[(j)\ is chosen in the form 

[0] = cxp { - ^ ^ [4>l (r) (4.9) 

+ (a/A)' [0,(T + A)-0,(T)]')rfr} 

We assume, as for the continuous limit, that (t) = 
(t + /3) ■ Of course, the representations (|4.6p and (|4.7I 
14. 9p are equivalent in the continuous limit. In Eq. (|4.8p 
one integrates the exponent over each slice instead of 
taking it in the middle of the slices as it is done in the 
second of Eqs. (j4.6|) . 

The function (/)r (r) is discontinuous and its typical 
form is represented in Fig. [51 The fluctuations of this 
function are essentially dependent on the ratio a/A. In 
the limit a ^ A, the typical average value {[<^>r{T)]'^)^r 



calculated with Wa, Eq. 



0r (t) 



, equals 



= Vo/A, 



(4.10) 



w 



whereas the correlation at neighboring slices is 



4>r (r) 4>r (r + A) 



Foa7A^<( <j)r{T) 



w 
(4-11) 

which shows that {[(l)r{T)~4>r{T+A)]'^)^ and ([0r(''')]')^^- 
are of the same order. 



In the opposite limit, a ^ A, the fluctuations of the 
difference of the fields at different slices are strongly sup- 
pressed and one can replace the difference in Eq. (|4.9p by 
the derivative. Then, we can use Eq. (|4.5p for estimates 
to obtain 



(r) 



Vo/ (2a) , 



(4.12) 



w 



while the average squared difference of (j) at neighboring 
slices equals 



; (r) - (r + A)) ) = A/a' « / (r)] ^ 

(4.13) 

Eqs. (|4?T2l |4T3| show that in the limit a > A the field 

(t) changes slowly from slice to slice. 

Now we are prepared to discuss the origin of the sign 
problem and how it can be avoided using the bosoniza- 
tion. 

Although the HS field is now discontinuous, Eq. (|4.8p 
is a standard continuous time representation for the par- 
tition function. This means that one can use Eqs. (|2.21I 
I2.25P replacing everywhere the function 0^ (r) by 0^- (''") ■ 
For example, the equation for the Green function takes 
the form 



-^-/^M)Gt"t(r,r') 



<5.,r'5(r-T'),(4.14) 



= — [l! ~ (7U(j)r (t) . 

Substituting the solution of Eq. (|4?T4l) into Eq. (|2?2T|) 
one comes to the fermionic determinant for Zf[(j)], 



det 



OT 



(4.15) 



In the same way as in Ref. [22 we can reduce the func- 
tional Zf[(j)] to the form 



ZM 



det 

r.a 



: det 

r.a 



l + Tr exp 



h[4>]dT 



N 



1 + Y\_cxp (- {ir - m' - f^0r-,;) A) 



1=1 



(4.16) 



In Eq. (|4.16p . the multipliers are ordered in I, such 
that the multiplier with Z = 1 is on the right. The 
same representation is usually used in the standard 
MC simulations [Actually, one replaces each multi- 
plier in the second line of Eq. (|4.16p by the product 
exp (— (e^ — /i') A) exp (cr0r,; A) , which is clearly justi- 
fied in the limit A ^ /?]. Eq. (|4.16p was obtained in 
Refs. [2^125! as a result of discretizing the time. At the 
same time, we see from Eqs. (j4.14[ I4.15P that only the 
HS field (pr {t) is discrete now, while the derivative d/dr 
is continuous as before. This allows us to bosonize the 
fermion model even in this discrete time representation. 
However, before making the proper transformations, let 
us discuss the origin of negative signs of Zf[(j)\. 
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C. The negative sign problem for fermions 

In the static case, when (jj^^i does not depend on /, the 
product in Eq. (|4.16p is triviahy calculated because the 
multipliers commute with each other. Then one comes 
to the standard expression for the partition function of 
fermions in a static external po 

tential, which is clearly positive. 

If the function 0^,; varies sufficiently slowly in time, the 
adiabatic approximation can be used and the function 
Zf[(l>] takes the form^^ 



N 

Zf[4>] ~l[[l + cxp(^-AY,^ 



(fe) 



(4.17) 



1=1 



where Ap'' is the eigenencrgy of A:-state of the Hamilto- 
nian obtained on each slice I. This expression is, 
again, positive. 

Actually, the negative sign of Zf[(j)\ may arise only 
when the function 0^,/ strongly fluctuates in time. 

Let us see how it happens. For a given (/)(r), we con- 
sider a matrix element Pnm of the product in Eq. (|4.16p 
for one of the spin directions 



p — 

^ ri m. 



where w)^™'' is the m-th eigenfunction of the Hamiltonian 
hi {4>) at the l-th slice [we assume that the eigcnfunctions 
D^™-* are real and the product in Eq. (|4.18p is ordered in 
time from the right to the left]. 

Using the completeness of the sets of the wave func- 
tions v''^!^ for each slice Z, 



E 



„(") 



N 



JJexp f-/i;A 



1=1 



(4.18) 



spins "up" and "down" [iVTj"*"' is the number of the eigcn- 
functions corresponding to the total number of the sites 
in the system]. The matrices P± are real but not neces- 
sarily symmetric. 

For fields </> (r) slowly varying in time, the eigenfunc- 

tions v). I on the neighboring slices are almost equal to 
each other and the matrices tt; are close to the unity ma- 
trix. In this case, the matrices tt; can be calculated using 
the standard perturbation theory for wave functions 
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where (5(/)r,; = (t>r,i+i — 4>r,i is small. Then, the off-diagonal 
elements of the matrices tt; take in the main order in 6(j) 
the following form 



2 i-^k-^Ti 



m ^ n 



l) 



(4.24) 
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In the limit A <C a, the characteristic values of the 
quantities |<50r,;| follow from Eqs. (|4. 12114. 13"]) 



\<f>r,i\-iVo/a) 



1/2 



(VoA/a 



2\l/2 



(4.25) 



while the average values of (f>r,i and (S(/>r,/ vanish. It is 
clear from Eq. (j4.24|) that |(50r,;| ^ \4'r,i\ and this justi- 
fies using the perturbation theory for computation of the 
matrix elements tt™", Eq. (|4.24p . 

As wc assume that ^ {tr,r',Vo}, the wave func- 
tions ai'c almost localized on the sites r„ and the 
states n are characterized by the positions of centers of 
the localization. Then, we estimate the denominators in 
Eq. (|4:24)) as 
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(fc) (fc) X 
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we represent the matrix element Pnm in the form 

N N 

Pnm = J2 5n,kjk,m CXp ( - A ^ [] ^f' + ^"' , 



{kl} 
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where 



fc, + ifei _ + (ki) 



(4.20) 



(4.21) 
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(4.26) 



Completely localized wave functions cannot give a contri- 
bution to non-diagonal matrix elements {54>r,i)™'^ in Eq. 
(|4.24p and one should take into account their overlap at 
finite try leading to finite values of the wave functions on 
sites in the neighborhood of the center of the localization 
tq. These values can be estimated [considering the tun- 
neling term try as a perturbation] as tr^ro \4>r,i — 4'ro,i\ ^ 
and one comes to typical absolute values of randomly 
fluctuating off-diagonal elements of the matrix tt/ 



.1/2 



m ^ n 



(4.27) 



are elements of orthogonal matrices tt;, and write the 
function Zf[(p] as 



det [l + P+]det [1 + P-] 



(4.22) 



where iV*°*'^' x iV^°*°' matrices P± are matrices with real 
matrix elements (P±)^„, Eq. (|4.18p . corresponding to 



The estimate (|4.27p is written under the assumption that 
all the states are non-degenerate, which is most proba- 
ble for randomly chosen fields (f>r,i- The matrix elements 
TT,""', Eq. (|4?27|) . are small in the limit A ^ but the 
matrix P, Eq. (|4.20p . contains a product of iV = /?/A 
matrices tt; and can essentially be different from a diago- 
nal one. Since the matrix elements tt; fluctuate randomly 
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with a gaussian distribution, we estimate the ofF-diagonal 
elements of the product of N matrices entering Eq. (|4.20p 

/ \ 1/2 

as ( r'/^o ) ■ This is a smaU number in the hmit of 
sufhciently high temperatures 



1/2 



(4.28) 



and Eq. (|4.28p determines the adiabatic hmit. In this 
hmit one comes to Eq. (|4.17p and the fermionic deter- 
minant Zf[(j)\ must be positive. AUernatively, one could 
take lower temperatures but large a ^ 13, which would 
lead to the same result. 

Let us discuss now what happens if a ^ /3 and the 
inequahty is not fulfilled. Using Eq. we 

write det P for one of the spin directions in the form 



AT N^"'"' 



N 



detP = cxp(-A^ Aj''Mj|det 



;7r;. 



(4.29) 
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As the matrices ni are orthogonal, one comes to values 
det TT; = 1 for all / and to positive values of det P in Eq. 
(|4.29p . Of course, this does not necessarily mean that 
the function Zf[(t>\, Eq. ()4.22|) . is positive, too. 

Matrix P given by Eq. ()4.20|) is real but not sym- 
metric. Such matrices cannot generally be diagonalizcd 
but one can always find a Jordan normal form for them. 
Assuming that P is an M x M real matrix one writes 

P = QJQ-\ (4.30) 

where the block diagonal matrix J can be written as 



J 



/Ji 
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and 



/ 1 \ 
p,; 1 ' 



1 

V 



The sizes Ki of the blocks Ji depend on the form of the 
matrix P but 



K, = M. 



For diagonalizable matrices each block in the matrix J 
is just a number and all Ki = I, such that n = Al. The 
numbers pi are solutions of the equation 



det {phi -P)^0, 



(4.31) 



where Im is the unity matrix, and are not necessarily 
real. However, as the matrix P is real, for any solution pi 
of Eq. (|4.3ip its complex conjugate p* is also a solution. 
One can see easily from Eqs. (|4.30[ I4.3ip that 



det[l + P]=llil+p,f\ detP = Y[i 



(4.32) 



As we have mentioned previously, in the adiabatic limit 
(|4.28p . when the function (j)r (r) is smooth in time, one 
obtains 



JV 

n 

1=1 



and one comes to inequalities pi > for all i. Deforming 
the function 0^ (t) , such that it can vary fast in the inter- 
val [0,/3], makes the matrix Jli=i essentially different 
from unity and some of the eigenvalues pi can become 
negative, which is a necessary condition for obtaining 
negative values of det [1 -f P]. 

In the limit A ^ a, all the matrix elements of the ma- 
trices TTi can be typically of order one and the changes 
of Pi can be large after a minimum change of a configu- 
ration of the field 0^,;- In this case nothing prevents pi 
from jumping from positive to large negative values and 
making dot [1 + P±\ negative. 

In the opposite limit a 3> A, a smooth deformation 
in time of the field results in a smooth motion of 
the eigenvalues pi. If we start with (f)rj slowly varying in 
the interval [0, (3], all pi are real and positive. Deforming 
the function 0^,; results in a motion of the parameters pi 
along the real axis. They cannot cross zero because this 
would mean det P"^-^ = 0, which is impossible because 
detTT; = 1 for all I. 

However, two different values pi and pk can collide and 
move to the complex plane forming a pair of complex con- 
jugated values p and p*. Again, this cannot change the 
sign because these two eigenvalues lead to a multiplier 
|1 in the determinant. A possibility to change the 

sign of the fermionic determinant Zf[(p\ first arises when 
the values p and p* collide again on the real axis at a point 
Po < 0. After this event one obtains again two real neg- 
ative eigenvalues pi and pk- Then, they move separately 
and one of them may cross the point —1 leading to a neg- 
ative value of det [1 -I- P±\. Still, this does not necessarily 
mean that the full fermionic determinant, Eq. (|4.16p . be- 
comes negative because the product of the contributions 
corresponding to spin "up" and spin "down" should be 
taken but the positivity of the fermionic determinant is 
no longer guaranteed. 

We see that, when deforming continuously the function 
0r,i, a chain of several events must happen before the 
fermionic determinant becomes negative. Therefore, the 
probability of negative values of Zf[(j)] may be reduced in 
the continuous limit, A <C a, with respect to the opposite 
one, A ^ a, but one should not expect that it vanishes. 
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Of course, using these arguments we cannot check 
whether the sign of Zflcj)] changes or not for a given 
function (jj^^i [unless it is sufficiently smooth in the entire 
interval [0, /3]] and this can be done only numerically. Ap- 
parently, negative values of Zf[(j)] are generally possible 
even in the limit A/a — > but we do not intend to clar- 
ify this question here. What we want to do is to replace 
the function Zf[(j)] by another function Zb[(j)\ which is 
positive by construction. In doing so we expect that, al- 
though these functions are generally different for a given 
4>r,i, replacing Zf[<j>] by Zb[(j)] in Eq. (|4.7p will not change 
the result of the integration over in the limit of a large 
number of sites in the system including both the main 
system and the bath. In the next subsection, we present 
arguments justifying this expectation. 



D. And how it can be avoided in the bosonic 
representation 

We start our discussion by generalizing the model un- 
der consideration to a model with an arbitrary complex 
electron-electron interaction writing instead of Eqs. (|2.41 



Z' = Tr exp -/SH (s) , H = Ha + Mnt («) , (4.33) 



where the operator Hq is given by Eq. (|2.2p with a 
properly shifted chemical potential and the interaction 



(4.34) 



replaces iJ^) in Eq. (H^. The variable s is an arbitrary 
complex variable and the original partition function Z, 
Eq. (f2^ . is obtained from Eq. (|4.33p by putting s = 1. 

For a finite system with a total number of sites Nd, the 
number of fermions cannot exceed 2Nd and the function 
Z'^ is a finite sum of exponential functions of s. There- 
fore, Z** is analytical everywhere in the complex plane of 
s except the infinite point. The function Z'^ decays in 
the left half-plane when Re s — > — oo but grows in the 
right half-plane when Re s — > oo. 

If the function Zf[(j>] in the integral in Eq. (|4.7p be- 
comes negative for sufficiently many configurations 0,-^;, 
making approximations in the integrand is dangerous be- 
cause the latter can be a fast oscillating functional of 
4>r,i- A simpler task is to calculate Z'' approximately for 
real negative s < because the fermionic determinant 
is strictly positive in this case. These values of s cor- 
respond to an attraction between the fermions and the 
sign problem does not exist in this case. As soon as the 
function is known for real s < one can try to make 
an analytical continuation to s = 1. 

For practical computations it is more convenient to de- 
couple the attraction term by a HS transformation cor- 
responding to an external "potential" rather than to a 



"magnetic field" acting of the spins used here because 
the latter enters the effective action with ^/s^ which is 
imaginary for real s < 0, and such a representation is 
not very convenient for numerics. However, we do not 
change the decoupling scheme for the present discussion 
and simply replace cf>r^i -> y/s4'r.i in Eq. (|4.16p . Then, 



the fermionic determinant Z 
negative s as 



can be written for real 



Zf[ 



dot 
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1=1 



(4.35) 

and it is clearly positive. 

Substituting Eq. (|435)) into Eq. (gj]) we obtain an 
integral with a positive integrand. Absence of oscillations 
of the integrand makes approximations more reliable and 
we apply the bosonization scheme for the real negative s 
simply substituting everywhere (j) by iy/—S(j). As a result, 
we write a new bosonic partition function Z^ [(jj] for real 
s < in the form 



Zo 



exp 



-y 

2i ^ 
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where the function A satisfies the equation 



(4.36) 



Hr^r' (t) Aa-r,r' {z) = -iu(7\^nry .cr^r,r';a (t) i?Q 

(4.37) 

with 



Hr^r' (t) ~ Ai (^ir ~ £r' — iuG \/~~S^ ^ y (t)^ 



-iA2— + A7 

OT 



(4.38) 



following from Eqs. ([2321 ElOl) . 

One can easily see after summation over a that the 
exponent in Eq. (|4.36p is real for real s < and the 
function Z^\(^ is positive. As the transformation from 
the fermions to bosons is almost exact and both the func- 
tions ^|[(^] and Zl\^ are positive, we expect that they 
should be close to each other for real s < 0. The in- 
troduction of the finite parameter 7 and of the bath is 
important for such s, too, because one should exclude 
"parasitic solutions" for Aa-r,r' (t). 

So, we have bosonized the fermionic system with at- 
traction for which the sign problem does not exist any- 
way. Integrating over (pr (t) in Eq. (|4.7p with Zj [(j>] and 

Z^[(j3\ we thus obtain two partition functions Zj and Z^ 
on the entire semi-axis of real s < and argue that they 
should be close to each other because Zj[(j>] and Z^[(j>] 

are. The functionals Zjlcj)] and Z^[(j>] are positive and 
a small difference between them should not result in a 
considerable difference between the functions Zi and Z^. 

In order to obtain the partition functions Zf and Zi, 
corresponding to the original model under consideration 
we perform an analytical continuation of the functions 
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Zj and in the complex plane of s from real negative 
s < to s = 1. In order to analytically continue Z^ one 
can use Eqs. ()4.331 14.34p giving a finite sum of exponen- 
tials of s and the analytical continuation is performed 
just putting s = 1 everywhere. Then, one obtains im- 
mediately all standard formulas for the fermionic system 
with repulsion and Zf[(j3\ derived from these formulas can 
be both positive and negative depending on 0^ (r) , which 
leads to the sign problem. 

Having constructed the function Z^ [cj)] for real negative 
s, Eqs. (|4. 36114. 35)) ■ one should integrate over (jj^- (r) using 
Eq. (|4.7|) and obtain the partition function Z^. After 
this function has been found for all real s < one can 
find uniquely an analytical continuation Z| for complex 
s coinciding with Z^ for real s < 0. This means that 
one can represent Z^ in a form of a convergent series in 
powers of the complex s in a certain region of complex s 
including the negative real half-axis. 

Below are some arguments for the analyticity of Z^. 
As we have seen, the difference between the terms of this 
expansion and those for the fermionic partition function 
Zj arises from a different treatment of two or more elec- 
tronic Green functions at coinciding Matsubara frequen- 
cies and momenta. This difference should be small in the 
limit of a large bath and small parameter 7 and therefore 
the Taylor expansion of the function Z§ should converge 
for all complex s in the same way as the expansion for 
Zj does. In Section IIIII we have checked explicitly the 
agreement between the partition functions Zj and Zh in 
the first two orders of the expansion in the interaction. 
The calculation can be repeated without changes for arbi- 
trary complex s including real s < 0. Thus, the function 
Z^ is analytic for all complex s except the infinite point 
and the physical partition function Zh can be obtained 
putting s = 1 in this scries. This result allows us to ex- 
pect that the function Zt is a good approximation for 
Zf- 

Having put s = 1 everywhere in the series for Z§ we can 
represent it again in a form of an expansion of the integral 
over (j) with [(/)], Eq. (|2.43p . thus coming back to the 
bosonization formulas obtained in the previous sections. 
This means, that replacing the original fermionic model 
with the repulsion by the bosonic one, we expect a precise 
agreement between the results obtained by integration 
over (f> of the functional Zf[(j)\ and Zii[4>] in the limit 7 — >■ 
and Nd — > 00. Although both the functionals can equally 
be used for analytical calculations, the use of Zhlfj}] may 
be preferable for the MC method because it is always 
positive in contrast to Zf[(j3\. 

The analyticity in the the complex plane is evident for 
the partition function Z", Eq. (|4.33p . and our arguments 
justifying the approximation done in the bosonization 
procedure have been based on this property. However, 
the analyticity of Z^ and Z^ does not automatically im- 
ply the analyticity of the functions Zj[(j)] and Z^[(f)] for 

any given (pr (r) and one should check this property us- 
ing explicit formulas. These functionals are close to each 



other for real negative s and would remain always close 
for other complex s including s = 1 if they were an- 
alytical in the complex plane of s for all (f>r (t), which 
would imply a convergent series in s. If this is not so one 
can expect that the functional does not well ap- 

proximate the functional [0] for certain cj) (t) in some 
regions of complex s including s = 1. In this case, the 
functional Zi\(j)\ used for our explicit calculations can be 
considerably different from the analytical continuation of 
from real s < to s = 1 and, hence, from Zj:[(j>]. 

In particular, the functionals Zh[0] and Zjicj)] may have 

different signs for certain functions 4>{t) when Zf[(j)\ be- 
comes negative. 

It is not easy to consider the analyticity properties 
of the functionals Zj[(j)] and Z^[(j)] analytically contin- 
ued from Eq. (|4.35p and Eqs. (|4. 36114. 35|) . respectively. 
It seems that Zjlcj)] is still analytical everywhere in the 
complex plane of s because the fermionic determinant 
is a finite product of exponentials and therefore can be 
expanded in series in integer powers of s. However, the 
same cannot be said about Z^ [(/>], Eqs. (|4. 36114. 35|) . and 
we think that this quantity cannot be analytical every- 
where in the complex plane of s for any 0^ (t). 

A simple way to see this is to try to repeat the regular- 
ization of Appendix El for complex s. The only difference 
with respect to s = 1 in the derivation is that the eigen- 
value in Eqs. (|2.47p written for complex Adr,r' is not just 
A-^*, where is is taken from Eq. (|2.46p . but generally 
another value A''^*. It coincides with A^* only for real 
Mry corresponding to positive real s. 

As a result, one should replace in Eqs. (jA3[ lASp the 
value A* by A*. Then, one cannot exclude that for some 
functions 0,- (t) and some complex s the eigenvalue E 
turns to zero. The subsequent integration over u in Eq. 
(|4.36p would lead to cuts in the complex plane of s. 

In this situation, the functional Zb[(f)] obtained by just 
putting s = 1 in Eqs. (j4. 36114. 35|) can differ from the 
analytical continuation from real s < to s = 1. Then, 
the functionals Zf[(f>] and Zt[4i] can essentially be differ- 
ent from each other for certain functions (f>r (t) and can 
even have different signs. At the same time, the physical 
partition functions Zf and Zh can still be close to each 
other because the functions Zj and Z^ are analytical ev- 
erywhere in the complex plane. 

This discussion may help the reader to understand why 
replacing the not necessarily positive fermionic determi- 
nant Z f[4>\ by the always positive functional Z},[(j)] we 
still expect the correct result for the physical partition 
function Z obtained after integration over all functions 
(j)r (r). 

The arguments of this subsection are strongly based on 
the assumption that the system is finite, which allowed 
us to speak about the analyticity of the partition func- 
tion Z^ in the interaction constant sV . This analyticity 
can be lost at low temperatures and real s < in the 
limit of an infinite sample when the system becomes a 
superconductor but this does not seem to create prob- 
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lems. A good quantitative agreement with exact resuhs 
is expected for large systems implying that they include 
the bath. 

Although the arguments presented in the last subsec- 
tions may help to understand the general scenario of what 
happens in the process of the bosonization, they are not 
rigorous and their final validity is to be checked numeri- 
cally. 

E. Some remarks about the regularization 

In this subsection, we want to emphasize that one can 
select the solution of Eq. p.32p in many different ways 
but only one of them, Eq. (|2.39[) . gives positive real Zi,[(j>] 
for all (f>r (t) . In principle, in the limit of a large number 
of sites in the bath all these solutions should lead to the 
same result and we simply choose the most convenient 
one. 

Of course, just taking the original electron Green 
function d^^f).^ (r, r') from Eq. (|2.24p we would satisfy 
Eq. (|2.32p for r' = r -|- i5 but we understand that using 
this solution, one encounters the sign problem. There 
are infinitely many other solutions of Eq. (|2.32p leading 
to non-positive weights in the MC procedure. 

To see this, let us carry out the regularization keeping 
small but finite (5 for r' = t + 5. Although the param- 
eter b is infinitesimally small, its presence can be very 
important for choosing the proper solution of Eq. (|2.32p . 
As we consider discontinuous functions 0^ (■''), taking the 
proper limit is especially important because t and r -t- 
can belong to different slices leading to completely differ- 
ent values 0r (■'') and (t -I- 0). So, we write Eq. (|2.26p 
in the form 

(|: - w (-)) G^:P:. ir,r + 5)=0 (4.39) 
with the operator M^'^} (z) , 

AfW (z) = 4 - eV - u<j^l+].,a (r) , (4.40) 

where 

¥+X,Ar)^4>rir)-^r'ir + 6) (4.41) 
and S +0. 

We can repeat the transformations of subsection III Bl 
keeping a finite 6 > and introducing the regularization 
with the parameter 7 and the bath. As a result, we come 
instead of Eqs. (|2.39[ 12.401) to the following equations 

-f^ly (^) i^) = -^<^riry,.¥+X,, (r) i?a, (4.42) 

where a = 1,2, and Hl'^J is a 2 x 2 matrix 

nl+] (r) = AiMI+J (z) + iA2^+ A7. (4.43) 



As we have seen previously, just putting 6 = and keep- 
ing the parameter 7 finite in Eqs. ()4.351 14.36|) results in 
strictly positive weights Zblfj}] because the solutions of 
the homogeneous equation in Eq. (|4.42p do not exist for 
finite 7 and one comes to a real solution for Aa;r,r' (z)- 
Putting 7 = for finite S we can choose the solution for 
Aa-r.r' {z) Corresponding to the initial fermionic problem 
and thus obtain the functional Zf[(j)\ that can be both 
positive and negative. 

In principle, one could obtain many other solutions 
of Eq. (|4.421 14.43^ keeping both 7 and 6 finite. Using 
the spectral expansion, Eq. (|2.52p . it is clear that the 
singularities in the function Ar^r' (z) can originate only 
from the zero eigenvalues = of the operator in 
the L.H.S. of Eq. (|2.32p . In order to avoid this problem 
we have introduced a regularization that discards such 
solutions. Actually, the crucial step in this regularization 
is taking equal times in the function $r.r' (■'') (putting 
,5 = in Eq. (lOTl) '). 

On the other hand, keeping a small (5 > makes the 
eigenvalues complex with an imaginary part of order 
d. Although this imaginary part is infinitesimally small 
in the limit (5 — > 0, it determines the imaginary part of 
the integral over u in Eq. (j2.43p arising near the poles 
corresponding to the zero eigenvalues E^ = 0. 

We calculate explicitly the imaginary part of the eigen- 
value E^ in Appendix [C] assuming that the eigenvalue 
E^ is very close to zero and show that it is finite for 
finite d but vanishes in the limit (5/7 0. 

Thus, our bosonization scheme should give different re- 
sults depending on how we treat Eqs. (j4.421 14.43p . The 
presence of the bath makes the partition function insen- 
sitive to the choice of the solution of this equation and we 
take the limit S/j ^ leading us to the the real positive 

Zb[4>]- 

In the next subsection, we discuss the basic proper- 
ties of Zi,[(t>\ and demonstrate that there is no reason to 
encounter the sign problem when using this functional. 

F. Properties of the bosonic action and its 
suitability for numerical investigations 

1. Bosonic action 

In this subsection we want to demonstrate that using 
Zf,[0] is convenient for computations and there are no 
singularities or bad features, like the negative sign, that 
would create problems. Final formulas will be brought to 
a form that can immediately be used for explicit compu- 
tations. We do not see any problems in taking the limit 
A — > 0, a/A — > in this representation and put a = 
in the beginning of the discussion. For convenience, we 
display in one place the formulas that may serve as the 
basis of computational schemes for the case of the on-site 
repulsion. 

We consider the system in a bath, which means that 
the HS field (f)r (t) is not equal to zero only in a region re- 
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stricted by a radius Rq (c.f. Eq. P-S?!) ). In the rest of the 
sample the particles do not interact with each other and 
the field 4'r (t) vanishes there. The interacting part of 
the free energy functional should be proportional to the 
number of the sites Nd located inside the sphere with 
the radius Rq. In order to obtain a good precision one 
needs to consider a large number of sites in the bath, such 
that the total number of the sites 7V*°*''' should consider- 
ably exceed Nj,. At the same time, considering the bath 
should not demand a large computation time because the 
HS field is identically zero in that region. 

We write the exact partition function Z of the original 
model in the form 



Z = lim 



/ ^b[0]WA[0]n;li#r,; 



N 
1 = 1 



(4.44) 



where (3 = iVA, 0^ (r) = 0^,; for - 1) A < r < ^A 
(see Fig. |6]), and Vq is the on-site repulsion. The values 
4>r^i satisfy the periodicity conditions 0^,; = 4'r.i+N- The 
weight Wa [0] has the form 



W^A [4>] = exp 



(4.45) 



in the region under the radius Rq. Outside this region 
the HS field is equal to zero identically. 

In order to write the correct functional Zi,[4>\ in the 
presence of the bath, the equations of Section |ll] have to 
be adjusted as concerns the shift of the chemical potential 
by the amount — Vb/2, Eq. p.l4p . We assume that main 
system is located in a potential well and that the height 
of the walls equals Vq/2. In other words, we shift the 
chemical potential /i everywhere in the space by Vo/2, 
although the electron-electron interaction is present only 
in the main system, using the value 

At' = M - VqI2 (4.46) 

This homogeneous shift is convenient because at /i' = 
the presence of the bath does not violate the particle-hole 
symmetry. Then, the functional Zfj[(j)\ equals 



Z}\(j)\ = Zoexp 



r.(T,a 



where 



and 



Zo = l[[l + cxp{ - I3ek 

k.a 



efc = (k) - + Vq/2. 



(4.47) 
(4.48) 

(4.49) 



In Eq. (|4.49l) . ji is the chemical potential of the original 
fermion model, Eq. (|2.2p . The functions A^.^ ^ (z) should 
be found from the linear equations 



7 Mr.r' 4 

Mr,r' - |: ^7 

-UUr^r'O'^r.r' (t) Ba, 



d_ 

dr 



where a = 1 , 2 and 



A'l , (z) 



B. = 



The operator M.r.r' and the function $r,r' (''') are defined 
as 



Mr.r' {z) ~ £r — £r' — UO'$r,r'(T), 
^r,r' (t) = ^r (t) - 0,-' (t) , 

Finally, the Fermi distribution 



E 



ik{r — r ) 



(4.51) 



(4.52) 



nk 



[exp{/3ea + 



Both the functions (pr {t) and Ar.r' {z)^ z ~ (r, cr, u), 
should be periodic in time: 0r (''') = 0r ("T" + P) and 



Aa:r,r' (t, CT, u) Aa-r^ (t + CT, u) 



(4.53) 



(4.50) 



The real parameter 7 should be taken as small as neces- 
sary to be sure that the final result is insensitive to its 
value. 

In order to calculate the partition function Z. one 
should solve Eq. (|4. 50114. 53)) for a set of functions 0^,; 
defined on the slices and a small parameter 7, substi- 
tute the solution AJ,.,. (z/) into Eq. (|4.47p and find the 

function Z\y\(^. This function is positive because the so- 
lution Aa-^r,r'{z{) must be real. Then, one can calculate 
Z in Eq. (|4.44p using as the probability in the MC 

sampling for finite A diminishing this length as much 
as necessary. We emphasize here that there should not 
be any principle limitation of the computational accuracy 
provided a sufficiently large system is considered and suf- 
ficiently small length of the time slices is used. 

These arguments are based on the fact that the solu- 
tions Ary iz) are real and non-singular, such that inte- 
grating over u does not generate imaginary parts. This 
property of the solution follows from the spectral expan- 
sion, Eqs. (|2.521 I2.54p supplemented by the proof that 
the real eigenvalues E entering the sum in Eq. (|2.52p in 
pairs E and ~E cannot be equal to zero (see Appendix 
E]) . Moreover, due to the fact that the eigenstates en- 
ter in pairs with the eigenvalues E and — _B, the entire 
exponent in Eqs. (|2.54[ I4.47[) does not become singular 
even for very small E due to the mutual compensation 
of the divergencies in the pairs. This follows also from 
the symmetry property of the Green function, Eq. (|2.53p . 
As all the parameters of the equations are real, there is 
no danger to obtain an imaginary part in the exponent 
in Eq. (|2.54p and therefore the weight is always 

positive. 



24 



In the limit of a very weak interaction Vq , typical (j)r (r) 
are small and one can neglect them in the expression for 
Mr,r', Eq. (|4.51|) . Then, eigenvalues Xk and eigenfunc- 
tions of this operator do not depend on the number 
of the slices and can easily be obtained in the form 



Ck (cos kr, sin kr) , = Ek, 



(4.54) 



where Ck are normalization coefficients. 

In this limit, one can solve Eq. (|4.50p by Fourier trans- 
forming both the sides in space and time. This allows 
one to write the Fourier transformed solution as 



K (k,k',c.) = 



au (£k - £k' - ('^k - n-w) '/)k-k', 
+ (£k - £k')^ + 7^ 



(4.55) 

Substituting Eq. into Eqs. I07)) and taking 

the limit 7 — > and A — > we come to the RPA result, 
Eq. p.l4p . The integral over u and summation over a 
are trivial in this calculation. 

The presence of the regularizing parameter 7 results in 
the vanishing of the contribution of the state with k = k', 
w = 0. However, the error is small because the sum 
extends over many states with different k that inevitably 
exist due to the bath. 



2. Integral form of the equation 



The differential with respect to time equation ()4.50|) 
can be written in the integral form with the help of the 
bare Green function G'^. ^ introduced as the 

solution of the equation 



5r,riSr' y^S (t - Ti) (4.56) 



The function ^, (''":''"i) is a 2 x 2 matrix and is 

uniquely defined by the boundary conditions 

^°r';ri,r; ('^''^l) = ^r,r' ;ri .rj ('^ + '^l ) 

= Gry;r,,riir,T,+/3) (4.57) 

The corresponding homogeneous equation does not have 
solutions due to the presence of the parameter 7. 

The bare Green function Q^^,.^^ ^, (t, ti) can easily be 
found by Fourier transformation and can be written ex- 
plicitly as 



Q° , ■ (r,Ti) 



(4.58) 



rp gik(r-ri)-ik'(r'-r'i)-iw(T-Ti) 

tj^total\^'^ wA2 4-(£k-£k')Al-h7A 
\ a I k,k 

where the Pauli matrices A, Ai, A2 are specified in 

Eq. ^rm . 



Eq. (|4.58|) is written for the entire system including the 
bath and the momenta k, k' correspond to this enlarged 
system. Accordingly, the total number of sites A^'Q*"' 
enters the normalization factor in Eq. (|4.58p . The state 
with a; = 0, Ek = £k' is included in the sum in Eq. (j4.58p . 

One can perform summation over the bosonic frequen- 
cies Lo = 27rTn in Eq. (I4.58|) using the formula 



exp {—iio (r — ri)) sgn (t — Ti) exp (—a (r — ri)) 
—iio + a I — exp I 



(4.59) 

where the symbol "sgn" stands for the sign and a is 
a number. This gives the following expression for the 
Green function G'^j.i.^^ ^1 (t, ti) for arbitrary coordinates 
and times 



-1 



exp 


6(t 


~ n) /n.rj 


2 sinh 







1 



= (£,-£>) A- 7A1 (4.60) 



where 



&(T-ri) = s.gn(r-Ti)/?/2- (r-n), (4.61) 



In Eq. (|4.60p the operators e^i and e^; act on the corre- 
sponding variables entering the J-functions on the right. 
The function g^y.^^y^ (r, n), Eq. (|4J0l) . is exphcitly real 
and free of singularities. 

Using the Green function ^, (t,ti) , Eq. ((TOI)) . 

we reduce Eq. ()4.50p to the form 

r 

rE / (^:n)u<T$,,,,;(ri)Ai 



I'l .r\ 



+ 5 (r - Ti) (5r,ri<5r'r; ^a;r^y^ (n , CT, u) cLti 
E / ^°r';ri,ri('r>n)Mnri,rlCrl>ri,rl {ti) BadTi 

(4.62) 



ri,r\ 



Eq. (|4.62p can serve as the basis of many numerical 
approaches. One can, e.g., solve this equation recur- 
sively neglecting in the zeroth approximation the func- 
tion ^ri,r{ {ti) in the L.H.S. and then calculating all or- 
ders of <&ri.rj (■'"i) or one can just invert the operator act- 
ing on Aa:r,r' (z) in the L.H.S. of the equation. Due to 
the rcgularization, the zero eigenvalues of the operator in 
the L.H.S. are discarded. Therefore all these operations 
are well defined and one should not encounter any singu- 
larity. In particular, one can invert the operator acting 
on Aa-r,r' in the L.H.S. of Eq. (gH]). 

As a result, after integration over u the exponent in 
the function Zh[(j)\ remains real and Zblcj)] positive. This 
demonstrates that this bosonization method of compu- 
tations is free of the sign problem. Substituting the so- 
lution obtained by inversion of the L.H.S. operator into 
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Eq. (|4.47p , the function [0] is given after discretization by 
I 



l,h 



du 



{r,r},{ri,r{} 



(4.63) 



r 



In Eq. ()4.63p the symbol Tr stands for the trace of the 
2x2 matrices introduced in the regularization, Eq. (|2.40p . 
The expression entering the square brackets under Tr is 
a 2 X 2 matrix in this "regularization space" and an N"^ x 
N"^ X N X N matrix with the space {r, r'} , {ri,r[} and 
time r, Ti matrix elements. This matrix is obtained from 
the matrix Q^a^ Ai using the conventional rules of matrix 
operations. The only peculiarity is that the coordinates 
enter in pairs {r, r'} and {ri,r'i} . Matrix elements can 
be written explicitly as 



(4.64) 

Eq. (|4.63p can directly be used for computations. The 
integral over u cannot lead to singularities, which is guar- 
anteed by the absence of the zero eigenvalues of the op- 
erator in the L.H.S. of Eq. (|4.62p written for an arbitrary 
u. This property is valid for an arbitrarily large <f>. 

In principle, one could easily perform the integration 
over u already analytically before doing numerics, which 
gives in the continuous time limit 



Zb[4>] = Zoexp 



X Tr 



(g"l>Ai)-ilnfl- (g"l>Ai)2 



1 rl3 rP 

O Y / f^T^n^r (t) nri.r; 

Jo Ja 



T.Ti 

{r,r},{ri,rj} 



(4.65) 



The integrand in Eq. (|4.65p is free of singularities as well. 
At small (j)r (t) one can expand the logarithm and take 
the first order of the expansion. This leads us again to 
the RPA formula, Eq. (|XTi|) . However, Eq. looks 
more convenient for MC computations because it allows 
one to peform updating more easily. 

Note that in the present paper, we decoupled the in- 
teraction by the Gaussian integration over (pr {t). For 
the on-site interaction, one can also use the "Ising spin" 
auxiliary field Sr^r^- All the steps of the derivation of 
the functional Zi)[(j)] can be repeated with this field. The 
final expression, Eq. (|4.44p . should be replaced by the 
following sum over all configurations of s^,; 



(4.66) 



where TV^j is the number of the sites in the system. The 
functional Zb [0] for the Ising spin field has the same form 



as previously with the field (pr (t) = Asr,;/A, coshA = 
exp(yoA/2), for (/ - 1) A < r < ^A. Using the "Ising 
spin" auxiliary field can further improve the efficiency of 
the computational schemes. 



3. A numerical test 

As we have discussed previously the final formulas, 
Eqs. (|i^liT7)) or Eq. are different from the for- 

mula of the original fermionic model, Eqs. (|4. 13114.1"^ . 
These deviations are due to the introduced bath and 
states effectively discarded by introducing the parame- 
ter 7. So, generally speaking there is no reason to expect 
that the functions Zf[fp] and are close to each other. 
As we have argued, they should be considerably different 
for, e.g., piece- wise functions with big jumps. The agree- 
ment between the partition functions Z is possible only 
after integration over (f)r and only in the limit of a large 
number of the sites iV^"*"' in the system. 

However, what one can expect is the agreement be- 
tween the values of Zf[cj)] and Zi,[(j>] for static 4> (one slice 
in time) for large Nd- Although we have demonstrated 
this agreement analytically in Appendix [B] it is instruc- 
tive to check this property numerically using Eqs. (|4.62[ 
I4.47P or Eq. (|4.63p . Such a test can also give a feeling of 
how well the numerical procedure based on using these 
equations can work. 

We present in Fig. [7] an example of the convergence be- 
tween Zf and Zb in a bath for the static case of one time 
slice. We consider a two sites Hubbard model embed- 
ded in a bath of Njf*"-^ — 2 sites with periodic boundary 
conditions. The HS fields (jjr are not equal to zero on 
the sites r = {1,2} but are equal to zero elsewhere, such 
that only the sites 1 and 2 are interacting. Both analyt- 
ical and numerical integration over u in Eq. (j4.63p are 
suitable and we use both of them. Integrating over u 
in Eq. (|4.63p analytically makes the computation faster 
but wc have found a good agreement between the results 
obtained using both the approaches. 

In Fig. [7] (a) the convergence between log (Zf/Zo) (up- 
per curve) and log {Zb/Zo) (lower curve) is demonstrated 
in the static the number of the sites in the bath is 

increased. Although the difference between these quanti- 
ties is fairly small, the complete convergence is not very 
fast. Therefore, in order to understand the dependence 
of the deviation log (Zf/Zb) on the size of the system we 
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FIG. 7: (a) Numerical evaluation of log(Z//Zo) (upper curve) 
and log{Zb/Zo) (lower curve) in the case of two sites Hubbard 
model imbedded in a bath of 7V^°*°' _ 2 sites in the static case 
(one time slice). The hopping parameter is taken as t = 4.0, 
the inverse temperature /3 = 1.0, the chemical potential fi' = 
0.1, and the regularization parameter 7 = 10~^. The HS fields 
on sites 1 and 2 are respectively (jjr^i ~ —1.0 and <j)r=2 ~ 
1.0. The entire system is a one dimensional ring in which the 
field (f) can be finite only on two neighboring sites, (b) The 
deviation \og{Zf /Zb) as a function of the inverse number of 
the sites in the system for the same parameters showing the 
1/7V^°'"' rate of convergence. 



plot in Fig. [7] (b) this quantity as a function of the in- 
verse total number of the sites Njj"*"-'- in the system. As 
previously, the HS is non-zero on two sites only and all 
the parameters are the same. 

Fig. [7] (b) shows that the deviation log(Z//Zt) de- 
cays with the size of the sample inversely proportional to 
the total number of the sites _/V*°*'^'. This dependence is 
natural because the regularization with the parameter 7 
cuts out iV^°*<»' of I^Nl°*-°-^Y states in the system. This ex- 
ample of the convergence completes the analytical proof 
of Appendix [Bl and shows that the bath is necessary in 
order to equate Zf and Z^. It shows as well that for any 
finite size system, Zf and Z^ are different, Z^ converging 
asymptotically towards Zf within an infinite bath. 

The numerical procedure carried out here can be ex- 
tended in a straightforward way to the case of the HS 
field (pr (t) varying in time, which corresponds to the 
case of many slices N . According to our analytical con- 
sideration there should not be principal problems with 
this general case and the function must remain real 
and positive. This is demonstrated in Fig. [S] 

For the calculation of the bosonic function the 
integral over u in Eq. (|4.63l) has been calculated analyt- 




FIG. 8: Distribution of (a) log(Z//Zo) and (b) \og{Zb/Zo) 
for 500 configurations for 4>r{T) randomly chosen according to 
the Gaussian distribution, Eq. (|2.12l) . for a system containing 
two interacting and one bath site in a ring geometry with 
parameters t — 1.0, /i — 2.0, Vo ~ 3.0, /3 = 14, and 7 = 
10~^ and for three time slices. Only log{Zf /Zq) can have the 
imaginary part leading to the negative sign of the fermionic 
determinant. 



ically (see Eq. (j4.65p ). The resulting matrix entering 
log{Zh[(j>]/ZQ) was evaluated using numerical diagonal- 
ization. This quantity remains real to a high precision 
[Imaginary parts due to numerical inaccuracy are of or- 
der 10~^^] showing the positivity of Zb[(j>], whereas the 
imaginary part of logZy[0] jumps frequently between 
and TT showing the sign fluctuations of Zf[<j)\. It is inter- 
esting to note that log ^/ is spread over a broader region 
of the real axis than log Zf, [0] . 

The form of the distribution of log Ztlcf)], Fig. |S1 shows 
that the sign problem is avoided in the bosonized rep- 
resentation. Moreover, the distribution is narrower than 
that of logZ/[0] and we hope that the computational 
scheme may be efficient. 



V. DISCUSSION 

We have mapped models of interacting fermions onto 
models describing collective bosonic excitations. This 
was performed by decoupling the interaction between the 
fermions by a Hubbard-Stratonovich transformation and 
considering non-interacting fermions in the fiuctuating 
HS field. The standard step of tracing out the fermions, 
which led to the fermionic determinant, was followed by 
converting this determinant into an integral containing in 
the integrand the solution of a linear differential equation. 
The latter step resembles writing quasiclassical equations 
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of Ref. [Ill but is now exact in the thermodynamic hmit, 
which has not been anticipated previously. 

The basic equation is actually an analogue of the von 
Neumann equation for the density matrix and our trans- 
formation can be interpreted by analogy with the replace- 
ment of the description of quantum mechanics in terms 
of wave functions by the description with the density ma- 
trix. The solution of the equation Ar^r' (t) plays the role 
of the density matrix. It obeys the bosonic boundary 
conditions Ar^r' (t) = ^r,r' (t + /?) and is a bosonic field. 

In order to regularize the obtained equations we as- 
sumed that the model of the interacting electrons is 
imbedded in a bath. The latter is a part of the sample 
where the electrons do not interact with each other. In 
other words, the sample we consider is in a contact with 
metallic leads. The introduction of the environment can 
serve as an additional support of the the analogy of our 
field Ar^r, (t) with the density matrix. 

We have checked the suitability of the new bosonic 
model for both analytical and numerical computations. 

As concerns the analytical approach, we have ex- 
pressed the solution of the equation for Ary (r) in terms 
of a functional integral over superfields ^' using the 
well-known in field theory Becchi-Rouct-Stora-Tyutin 
transformation^SiH. This allowed us to integrate out the 
HS field (j) and reduce the original fermionic model to 
a model of interacting bosons. One can work with this 
model using expansions in the interaction and we have 
compared the results of the perturbation theory with 
those for the initial fermionic model up to the second 
order in the interaction. This comparison demonstrated 
the exact agreement. 

In a short future, a renormalization group scheme anal- 
ogous to the one suggested previously in the quasiclassi- 
cal approacbi^ will be developed and applied for studying 
anomalous contributions to thcrmodynamical quantities. 
More complicated models for strongly correlated systems 
can also be studied analytically using the bosonization 
approach developed here. 

As concerns using the developed scheme for the MC 
calculations, the superfield theory obtained after aver- 
aging over (p with the help of the superfields ^' is not 
directly suitable for this purpose. Since we do not know 
how to reach the final form without introducing inte- 
grals over Grassmann variables, the only possibility for 
numerics is to calculate using the bosonic model in the 
fluctuating HS field. There is a very strong motivation 
to try to study the bosonic model using the MC method. 
It is well-known that using the MC method for fermion 
models suffers generically from the famous negative sign 
problem. Since we have mapped the fermionic model 
onto a bosonic one, one might have a hope that the sign 
problem will not show up. 

We have analyzed this idea and come to the conclu- 
sion that the MC simulations for the bosonized model 
shall be free of the sign problem. Practically, one can 
simply forget how the bosonic model has been obtained 
from the initial fermionic one and start investigating it by 



discrctizing time. Then one can make a HS transforma- 
tion with a discrete field (j)r (r) and obtain the formulas 
of subsection IIVFI This would be logically sufiicient for 
justifying our scheme. 

However, we have considered this issue in great detail 
in subsections IIV Bl and IIV CI When deriving the basic 
equation for the bosonic field Ar^r' {t) , the fact that only 
the first order derivative with respect to r enters the 
equations for the electron Green functions, Eqs. (|2.241 
I2.25p . is crucial. Numerical schemes are based on dis- 
crctizing the time but just replacing the derivative by a 
finite difference would make the reduction to the bosonic 
field Ary (r) impossible. We suggested a procedure of 
time discretization that enabled us to use a piece-wise 
HS field (f)r,i keeping the time continuous and carry our 
the bosonization. A crucial step in the transformation of 
the fcrmions into the bosons is the regularization based 
on the introduction of the bath and of the small param- 
eter 7. 

The importance of this regularization can be under- 
stood considering the conventional diagrammatic expan- 
sion for fcrmions. In this language, the bosonization 
means converting pairs of the fermion Green functions 
into a propagator of bosonic excitations. However, it 
is not possible to do this transformation if two or more 
Green functions have coinciding Matsubara frequencies 
and momenta. 

The regularization based on introducing a small pa- 
rameter 7 removes such states and the presence of the 
bath makes their contribution small in the limit of a large 
number of sites in the bath. Our mapping can be used 
for any temperature, interaction and dimension of the 
system and, in this sense, the mapping is exact. At the 
same time, the presence of a sufficiently large bath is nec- 
essary to provide the agreement between the boson and 
fermion models. 

We have discussed in subsections IIV CI and IIVDI the 
origin of the negative sign in the fermionic determinant 
Z f[4>\ and demonstrated that the corresponding func- 
tional Z}j[(f>\, being generally different from Zf[(j)\, is al- 
ways positive. We have done this representing the so- 
lution for the differential equation for the bosonic field 
Ar.r' (t) in terms of a sum containing cigenfunctions and 
eigenvalues, Eqs. p.52| I2.54|) . With the regularization 
scheme developed in Section [Hi all the quantities enter- 
ing these equations are real. In the language of Eqs. 
(|2.52l I2.54P , one would encounter the sign problem if an 
eigenvalue could turn to zero at certain u. Then, 
the integration over u near such poles might generate an 
imaginary part and make the function not neces- 

sarily real and positive. However, we have demonstrated 
(see Appendix|A| that the eigenvalues cannot turn to 
zero and there are no singularities in this spectral expan- 
sion. This allowed us to demonstrate that the distribu- 
tion Zi,[(j>] that should be used for MC sampling is always 
real and positive, hence showing that there is a chance to 
overcome the negative sign problem by the bosonization. 

Although the functionals Zflc/)] and ^^[0] may even 
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have opposite signs for certain functions 0, wc present ar- 
guments that after integration over (jj the partition func- 
tion Zb obtained in the bosonic model has to be a good 
approximation to the partition function Zf of the orig- 
inal fermion model provided both the models are taken 
with the bath and the regularizing parameter 7 is small. 
The absence of the sign oscillations in Zi,[(j>] should lead 
to smaller variations of the modulus of this quantity and, 
hence, to a better convergence in the MC simulations. 

The logics of all these manipulations can be illustrated 
calculating two different integrals 

/>oo />oo 

Ii = V2 cosax^dx, I e-'^'^^dx (5.1) 

Jo Jo 

Calculating the integral Ii analytically, one can reduce it 
by turning the contours of the integration to the integral 
I2 and show exactly that Ii = I2 for any a. 

At the same time, we can discretize the variable x, 
replace the integrals by sums and calculate these sums 
by the MC method. Then, one can see that difficulties 
in the calculation of these sums are quite different. The 
integrand of I2 is positive and the sum converges very 
fast. In contrast, the integrand of Ii oscillates faster 
and faster with growing x and one cannot use it as the 
weight in the MC sampling. Of course, one could replace 
the integrand by its modulus and use the latter as the 
weight but then one would calculate the average sign of 
cos ax^ as in Eq. (|4.3p , which is not a very pleasant task. 

Our scheme of the replacement of the fermionic deter- 
minant Zf[(j)\ by the bosonic functional Zh[(j)\ resembles 
the analytical reduction of the integral Ii to the integral 
12- This transformation can easily be carried out in the 
continuous limit by using the powerful theory of complex 
variables but making the same for finite length of slices is 
difficult and cannot be done exactly. It is clear that the 
discrete versions of the integrals Ii and I2 are different 
and the only justification for a replacement of Ii by I2 is 
that they must be equal to each other in the continuous 
limit. 

A nice example relating the occurrence of the sign 
problem to the representation frame of a physical prob- 
lem is discussed on p. 103 of Ref. In this example, 
the model of a single spin in an external magnetic field 
is discussed. Although the corresponding Hamiltonian 
H = — h • <T can be solved in a simple way, tackling 
the problem using the path integral approach can either 
lead to sums containing both positive and negative terms 
[choosing h = {hx,hy,0) to lie in the xy-place], or lead 
to purely positive matrix elements [after rotating h into 
the xz-plane, h = {hx,0, h.^)]. While the result of the 
former choice resembles the fermionic sign problem, this 
sign problem is apparently "solved" by switching to the 
other reference frame. 

This illuminating example shows that in some cases 
the sign problem can be overcome by a simple analytical 
transformation. Our approach for solving the fermionic 
sign problem looks similar: we make a rotation in a "gen- 
eralized space" of bosons and fcrmions from the "fermion 



plane" to the "boson plane" . Actually, it is not surpris- 
ing that one does not encounter the sign problem in a 
bosonic model, but the fact that there can be an exact 
mapping between a fermion and boson model has not 
been anticipated previously. 

The sign problem is a long standing problem of quan- 
tum MC computations for fermionic systems. There 
have been many attempts to overcome this problem 
writing sophisticated algorithms but it has not gener- 
ally been solved until now. It has even been asserted^ 
that the sign problem belonged to the class of NP-hard 
problems21, which would imply that chances to solve it 
are very low. However, we do not think that the argu- 
ments presented in Ref. [2^ can be considered as a rig- 
orous proof. Even if our bosonization scheme will really 
work for the fermion models, we do not know yet how to 
apply it for any NP problem. 

The functional Zb[(j)\ can be computed for any (j)r (r) 
using directly Eq. (|4.63|) . One can perform integration 
over u either numerically or analytically. The numerical 
integration over u looks more convenient because one can 
easily carry out the MC updating. Then, one can use the 
standard MC scheme for integration over the field 4>r (r) . 
One can also use the "Ising spin" auxiliary field^^ instead 
of the Gaussian field 0, which can make the computation 
faster. 

In order to make an independent check and get a feel- 
ing of how one could compute using Eq. (|4.63|) we have 
calculated the function Zi, [cp] for a static field (pr integrat- 
ing over u both analytically and numerically and found 
a good agreement with the corresponding function Zf[(j}] 
of the original fermionic model. We do not expect es- 
sential difficulties in extending the computation to time- 
dependent HS fields (pr [t) and hope that our bosoniza- 
tion scheme will be checked numerically in the nearest 
future. 



Acknowledgements 

Wc thank Transregio 12 of DFG, the French ANR for 
financial support and the Aspen Center for Physics where 
part of this work was completed. We acknowledge useful 
discussions with Y. Alhassid, A. Bulgac, D. Galanakis, 
M. Jarrell, S. Kettemann, M.Yu. Kharitonov, O. Parcol- 
let, M. Troyer, Ph. Werner and S-X. Yang. 

Appendix A: Absence of zero eigenvalues in the 
regularized model 

Here, we demonstrate that the operator 'Hr,r' {T)^ 
Eq. (|2.40|) . does not have zero eigenvalues . This 
property is guaranteed by the presence of the regular- 
izer 7. If we put 7 = in Eq. (j2.40p . zero eigenvalues E 
exist. 

Putting 7 = we come to equations ()2.46l I2.47P . As 
the operators in the L.H.S. are not hcrmitian, the eigen- 
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values are complex and can turn to zero for certain 
fields (f>r (r). 

Suppose that the modulus of an eigenvalue A becomes 
very small or can even turn to zero and let us show that, 
nevertheless, the corresponding eigenvalue E remains fi- 
nite for a finite 7. 

Let Wr,r' (''■) be the cigcnfunction corresponding to the 
eigenvalue A in the first equation in (|2.46|) . Another func- 
tion V* ^1 (t) should correspond to the eigenvalue A* in 
the second equation in (|2.47p . 

The small values of | A| correspond to small values of the 
eigenvalues E of Eq. (|2.48|) of a state with an eigenvector 
Sry ij). At 7 = and A = the solution of Eq. (|2.48p is 
very simple. One obtains immediately i? = 0, while the 
solution Sr.r' ij) can be written in the form 



Sry (t) = ci5,V (r)+C25,V (r) 
with the vectors 5° ^, (r) and S\ ^, (r) given by 



ir) = 







(Al) 



(A2) 



and arbitrary coefficients ci and C2 . This means that the 
state with £■ = is degenerate. 

The case of small A and 7 can be considered using 
quantum mechanical perturbation theory. As the state 
with £■ = is degenerate we seek for the solution Sr,r' (t) 
writing it in the form of Eq. (jAl[) . Following the standard 
scheme other states are neglected in this expansion. 

Substituting Eq. (|XT|) into Eq. (|2.48p we multiply both 
sides of the equation first by S^y (r) , then by S^*^, (r) , 
sum over r, r' and integrate over r. This gives us a system 
of two equations for the coefficients ci and C2, 



{E - 7) xici - C2A = 

-CiA* + {E + -/)X2C2 = 



(A3) 



where 



Xl 



X2 



V / \Vry {T)\^dT, 

V / \V„., (T)l^dT. 

/ Jo 



(A4) 



Non-trivial solutions of Eq. (|A3|) exist provided the de- 
terminant equals to zero. This condition gives us the 
eigenvalues E 



E = ±\h^ + 



lAI 



X1X2 



(A5) 



Eq. (|A5[) demonstrates explicitly that the eigenvalues 
cannot turn to zero even in the situation when A = 0. 
This is a well known effect of level repulsion. The ab- 
sence of the zero eigenvalues E justifies our method of 
the regularization. 



Appendix B: Checking the regularization for a static 
Hubbard-Stratonovich field 



In this Appendix we demonstrate that using the spec- 
tral expansion, Eqs. (j2.521 12.54p . we come in the case of 
static HS fields back to Eq. (|2.2ip . This can serve as a 
check of the transformations we have carried out in or- 
der to derive Eqs. (|2.521 12.54p and can help the reader to 
visualize our regularization scheme. 

Assuming that the HS field 1/)^ does not depend on time 
we write time independent solutions S^^, of Eq. (|2.48p in 
the form 



(Bl) 



CK 

'-'l;r.r' 
cK 
*-'2;r,r' 

where a and b are coefficients and the functions are 
the eigenfunctions of the operator hr, Eq. (|2.24p 

hrw^ = X'^w^ (B2) 

In the static case, only time independent solutions con- 
tribute into the exponent in Eq. (|2.54p and for such states 
K = {k,k'}. 

Substituting Eq. (|BT1) into Eq. (|2.48p we come to a 
system of two linear equations 



a - X" 
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E^'b 



(B3) 



= E'^a 



a7 + 6(A*^- A*^'') 
Solution of these equations leads to the eigenvalue 



(B4) 



Substituting Eq. (|2.52p into Eq. (|2.54p we have to calcu- 
late the sum over the eigenstates S*^^, and thus find the 
solutions for a and b from Eq. (jB3p using the eigenvalues 
E^, Eq. (|B4)) . The function Z [^] can be written as 



Z [(f>] = Zq exp 



IE 



E 

K 



( qK* oK , oK* qK \ 

\'^l-,r,r"-'2;ri,r[ ^ '-'2-r,r' ^ l-.n ,r[ J 



E^ 



udu 



(B5) 



The subscripts 1 and 2 in Eq. (jBSp correspond to two 
different components of the vector S*^^, , Eq. (|Bip . 

One can easily obtain from Eqs. (jB3[ IB4p that the 
product ab entering Eq. (|B5P takes the form 



ab = ±- 



I 



A*^ - A*^'' 



(B6) 



^y^(Afe-Afc')' + 7' 

and we come to the following form of the function Z [(jj] 



Z\ 



Zq exp 



-udu 



i(Afe_A'=')^+72 



(B7) 
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We make a further transformation of Eq. (jB7[) using 
Eq. (|2.26p and a corresponding equation written in the 
absence of the HS field 4>rcr ■ Subtracting these equations 
from each other we obtain 



(B8) 

we reduce the func- 



in this appendix are an extension of those obtained in 
appendix]^ to the case of finite 5. 

Assuming that 5 is smah we expand the field (f)r' (t + 6) 
entering Eqs. (|4. 33114. 35)) in this parameter and repre- 
sent the eigenvalue problem of the operator (t), 
Eq. (|4.36p . in the form 



(z) = ESr,r' {z) ■ 



(CI) 



Substituting Eq. into Eq. (|B3 

tion Z [(/)] to the following form 



(r) - Ur,r 



{z) + 8ry (2) , 



Z [(/)] = Zo exp 



E 



^0 



dudracf)'. 



where the hermitian operator "Hr-.r' (z) is defined in 
, fc'Eq. (|2.40p and the matrix dr,r' (z) can be written as 



(z) 



{T,r + 0)-Gi°l, (r,r + 0) 





OT 



auS 



d4>^, (t 





(C2) 



In the limit 7—^0 the state with k = k' drops out from 
the sum and we obtain finally using the orthogonality of 
the functions 



= Zq exp 



Jo 



0) 



T,T + 0) - 



r,ri ,r[ .a,k 



[g^:X„ (t, r + 0) - g(°;,, (t, r + 0) ) ] ] drd« (BIO) 



We see that, containing the second sum, Eq. (|B10[) 
is somewhat different from Eq. (|2.2ip . At this point 
we should recall that the system we consider contains 
a bath making the number of states large. The normal- 
izcd cigenfunctions are proportional to (7V^°*'»') 



Eq. (|Cip is a direct generalization of Eq. (j2.48p to the 
(B9)casc of the split times. We assume that the eigenvalue 
E is close to zero. The regularization with the finite 7 
leads to level repulsion. All eigenvalues calculated 
for 5 = remain finite and real. 

Considering finite 8 drastically changes the situation 
even if 5 is infinitcsimally small. This is a consequence 
of the fact that the matrix 5r.r' {z) is not hermitian and 
therefore the operator (r) acting on Sr,r' {z) in the 
L.H.S. of Eq. (jCip is not hermitian, too. This means that 
the eigenvalues E can now be complex. 

In order to calculate the eigenvalue E we use, as in ap- 
pendix|X]the degenerate perturbation theory considering 
the matrices 5ry {z) and 7A as a perturbation. 

We represent the solution of Eq. (|Cip in the form of 
Eq. jJH and multiply both sides of Eq. ((CT|) by 

5"°*/ (t) and then by 5*^*^, (r) , summing after that over 
r, r'and integrating over r. As a result, we obtain equa- 
tions for the coefficients ci and C2 



Therefore, the second sum in the exponent in Eq. (jBlOp 
is much smaller than the first term and can be neglected 
provided the total number of the sites in the system 
^totai -g jg^j.gg^ rpj^jg example of the static HS field 
demonstrates very well how the regularization works. 



Appendix C: Splitting times in the equation for 

In this appendix we demonstrate that keeping in 
Eqs. (|4. 33114. 36)) the parameter 5 finite one can obtain 
imaginary part in the eigenvalues E^ of the operator 
•H^-t? (t), Eq. gjS]). The finite 6 appears because the 
times r and r' entering the Green function g''^^} (t, r') 
should be slightly different, t' = t + 5. Although the 
value of 5 can be infinitcsimally small, the imaginary 
part of E^ can be very important in situations when 
the eigenvalue E^ is close to zero. The results derived 



^E - 7j xici 
- (A* - A2(51 



-C2 A + Ai(5 = 



ci + LB + 7 a;2C2 = 



(C3) 



with positive numbers xi and X2 given by Eq. ()A4p . The 
complex parameters Ai and A2 have the form 



Ai 



r.r' 



(r) 



dcj)r' (t) 
dT 

K (r) 



Vr,r' (j) dr, (C4) 



V* ^1 (r) dr 



where the functions w^.r' (t) and u^.r' (t) and the complex 
eigenvalues A are introduced in Eqs. (|2.461 1^115)) . 

In the limit |A| — > the functions Vr.r' (t) and Vr.r' (t) 
are related to each other as 



(t) 



Ar) 



(C5) 
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and we obtain 



A2 — Ai — A 



Substituting Eq. (|C6[) into Eqs. (|C3P and solving the lat- 
ter we obtain the eigenvalue E 



E ^± 



\ 



\\r + 25ilm XX 



T + 



X1X2 



We see from Eqs. (jC71IC8p that a finite value of 5 that de- 
termines the splitting of the times makes the eigenvalues 
(C6) E complex. The imaginary part, being infinitesimally 
small is very important for small A generating a finite 
imaginary part of the integral over u in Eq. (|2.54p . Work- 
ing in the limit of 7 ^ (5 one would inevitably encounter 
the sign problem. In contrast, putting (5 = and work- 
ing with small but finite 7 gives strictly positive Z [(j)] 
[Qf) a-nd one avoids the sign problem. 



Eq. (|C7|) demonstrates that the value of E near the 
point A = strongly depends on the ratio 7/^. In the 
limit |A| ^ (5 7 we obtain 



E = ± 



\X\+i5- 



Im l^XX 

w 



1 

(2:1X2) 



1/2 



(C8) 
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